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SYNOPSIS 


Viscous flow of Newtonian and non-Newtonian fluids past an 
assemblage of particles represents an idealisation of processes 
occurring in diverse engineering applications such as ground water 
hydrology, polymer processing, biotechnology and numerous particle 
separation techniques. Despite their overwhelming importance, the 
progress in the modelling of transfer processes in these systems 
has been rather slow. This is so partly due to the formidable 
difficulties encountered in describing the interparticle 
interactions satisfactorily and mathematical complexities arising 
due to non-linear governing equations The currently available 
theoretical studies concentrate mainly on evaluation of 
macroscopic quantities such as drag coefficients and Sherwood 
number. Detailed studies corresponding to intermediate Reynolds 
number flows are scarce for both Newtonian and non-Newtonian 
fluids. In fact, even for creeping flow of non-Newtonian fluids, 
the detailed flow field solutions have not been reported in open 
literature. 

In present work, comprehensive theoretical results on the 
detailed structure of the flow field, namely, pressure and shear 
stress distribution, angular and radial velocity profiles, 
streamline patterns as well as the macroscopic parameters such as 
drag coefficients have been obtained over a wide range of bed 
voidage, Reynolds number and rheological behaviour of the fluid 
The major contribution lies in providing physical insight for the 
observed drag behaviour based on the detailed flow field structure 
prevailing around the particle. 
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In Chapter I, the basic features of flow past a 
multi-particle assemblage are introduced and the scope of this 
study is outlined In Chapter II, the methodologies available for 
modelling interparticle interactions are reviewed and the 
pertinent literature is critically examined. Among the various 
approaches, the capillary and cell models are discussed at length. 
The serious limitations of the capillary model (apart from its 
strengths) are highlighted and the relative advantage of the cell 
model is brought forth The central concepts of the cell model 
regarding the division of the particle assemblage into several 
unit cells each housing a single particle and the subsequent 
analysis of flow through the assemblage in terms of flow within a 
typical cell are described. The widely accepted free-surface and 
zero-vorticity cell models are elaborately discussed 

In Chapter III, the application of free surface cell model to 
the steady flow of a Generalized Newtonian fluid past an 
assemblage of spherical particles is presented. The Generalized 
Newtonian fluid description provides a convenient representation 
of Newtonian and power-law type non-Newtonian fluids through a 
single rheological equation. The non-dimensional continuity and 
momentum equations in conjunction with appropriate boundary 
conditions have been solved by the FEM procedure based on the 
Galerkin's weighted residual approach. The discretization of the 
flow domain has been done using eight-noded quadrilateral 
elements. The nodal equations for velocity components and 
pressure have been derived respectively from the application of 
the tesidue principle to the momentum and continuity equations 


On the boundary, however, the nodal equations are obtained by the 
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direct implementation of the boundary conditions The resulting 
system of non-linear algebraic equations are solved iteratively 
using the FRONTAL t echn ique 

I*- Chapter I V , the validity and accuracy of the numerical 
scheme are examined. Comparison of present results with the 
analytical solution of Happel for the creeping flow of a Newtonian 
fluid shows a perfect agreement of drag coefficients (pressure, 
friction and total), surface stresses and the velocity and 
pressure field solutions. For intermediate Reynolds number flow 
of Nev'tonian fluid through a dense assemblage, the drag results of 
present study are in excellent agreement with the perturbation 
results of El-Kaissy and Homsy At higher voidage, the Reynolds 
number upto which agreement exists decreases For the creeping 
flow of a power law fluid, fair agreement is observed between the 
drag pi edictions o< present work and the upper and lower bounds 
available m the literature m the range of power law index, 1 £ n 
& 0.5. Experimentation with the radial step size and grid pattern 
indicates that the number of elements in radial direction required 
for accurate prediction increases with voidage; also, for dilute 
assemblages, a grid with exponential size variation is found to be 
far superior xo a uniform mesh Convergence of drag coefficients 
requires less iterations than those needed for the convergence of 
velocity and pressure solutions. The rate of convergence is 
reduced when the voidage and power law index are decreased, while 
the Reynolds number has no appreciable effect. 

In Chapter V, Newtonian flow results are presented over a 
wide range of voidage and Reynolds number. For dense 
assemblages , many creeping flow features such as inverse 
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proportionality between drag coefficients and Re, symmetric 
surface shear stress profiles, etc. are displayed even at 
intermediate Reynolds number. The critical Reynolds number upto 
which the creeping flow solution is valid depends on the voidage 
(<£ ) . At low voidages, pressure drag dominates over friction drag 
(Cpp > Cpp) and the reverse is true at high voidages The switch 
over between C^p and C^p takes place around e = 0.7. Comparison 
of present results with the experimental correlation of Ergun for 
packed beds shows that the free surface cell model is accurate for 
Re < 10 and for higher Reynolds number (Re), the zero vorticity 
model predictions are slightly better, although both models 
diverge from Ergun’s correlation as Re increases For 
fluidization best correspondence between the predictions of widely 
accepted Richardson and Zaki correlation and the present results 
is observed at Re = 20.0. 

In Chapter VI, the drag and the detailed flow results for a 
power law fluid are presented. For creeping flow, the pressure 
and friction drag coefficients vary as a power function of the 
average Newtonian surface shear rate for the range of power law 
index, 1 > n > 0.5. With respect to voidage, the pressure drag 
dominates over the friction drag for dense assemblages and vice 
versa for dilute assemblages. For 0.3 < <e < 0.9, the ratio 
C,.^ /C_„ decreases as n is decreased, while an exactly opposite 
trend is observed for voidage close to single sphere limit. The 
surface shear stress and pressure profiles become increasingly 
flat as the degree of non-Newtonian behaviour increases; this 
trend is more evident at low voidage. The flattening of the 
surface shear stress is caused by the simultaneous decrease in the 
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fluid viscosity and the increase in the shear rate adjacent to the 
particle The variation of surface viscosity with angle is seen 
to be prominent only near the axis The effects of power law 
index and voidage at intermediate Reynolds number flows are 
essentially similar to those of creeping flow. The inertial 
effects on drag are noticeable only for high Re, high voidage and 
low values of n. The surface shear stress and pressure profiles 
are, however, significantly affected, shear stress becomes 
increasingly asymmetric about the mid-plane and the pressure 
approaches the potential flow pressure profile. Finally, present 
theoretical predictions have been substantiated by performing 
detailed comparisons with appropriate experimental data available 
in the literature for packed beds, fluidized beds and single 
particle falling through a fluid under the action of gravity 
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CHAPTER I 


INTRODUCTION 

1.1 MULTI -PARTICLE SYSTEMS AND THEIR APPLICATIONS 

Fluid flow through an assemblage of particles has received 
wide research attention both theoretically as well as 
experimentally due to its importance in several engineering 
disciplines. The flow configurations involving multi-particle 
assemblages can be classified In many ways based on factors such 
as the proportion of solid to fluid and the relative motion 
between the solid and fluid phase, etc. Depending upon the 
relative proportion of solid to fluid, these are defined as 
concentrated or dilute particle assemblages. In practical 
applications, three broad classifications can be identified 
according to the direction of motion of the particle assemblage in 
relation to that of the fluid flow, namely, sedimentation, packed 
bed and fluidization. In certain other applications where the 
specific surface of the solid matrix is very high and pores 
comprising the void space are relatively narrow, the system may be 
termed as a porous medium. Each of these configurations of 
multi-particle assemblages Is of considerable practical 
signl f icance 

The flow of fluids through packed beds is an operation 
encountered in various engineering applications. Common examples 
concerning single phase flow are: 
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(a) Filtration of industrial suspensions , water treatment through 
sand filters, etc. 

(b) Mixing by means of flow through a porous bed in order to 
equalize temperature, velocity field or composition in the 
fluid. 

(c) Flow through a porous bed in underground hydraulics . 

The process of sedimentation is commonly employed in chemical 
and petroleum Industries as a way of separating particles from 
fluid or separating particles of different sizes. Examples of 
such separations include dewatering of coal slurries, 
clarification of waste water, removal of sugar crystals from syrup 
and the processing of drilling and mining fluids containing rocks 
and mineral particles of different sizes. The separation of 
different particles by sedimentation is also the basis of some 
laboratory techniques for determining the distribution of particle 
sizes in particulate dispersion. 

The process of fluidization of solid particles by a liquid 
which was considered to be just a scientific curiosity, is rapidly 
becoming more and more important from industrial point of view. 
New processes involving liquid-fluidized beds are being developed 
in areas such as hydrometallurgy, food processing, biochemical 
technology, water treatment, etc. The operations Involved may 
include crystallization, 1 1x1 vat ion , ion exchange, adsorption, 
enzyme catalysis or cell culturing . 

Innumerable efforts have been made to analyse the Newtonian 
flow through multi-particle assemblages and detailed reviews are 
available in several monographs. The classic treatise by Happel 
and Brenner (1965) outlines a rational approach to analyse the 
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fluid particle dynamics based on principles of fluid mechanics. 

In the last three decades, wide interest has been evinced in 
the study of non-Newtonian fluid flow in multi-particle systems, 
as it constitutes an important part of many industrial processes. 
Typical examples Include flow of oil through porous rocks, polymer 
flooding in tertiary oil recovery operations, filtration of 
polymer solutions and slurries, ore-processing, catalytic 
polymerization in hydroxydation process, waste water treatment and 
ion exchange beds which are used to recover metals like uranium 
from slurries and sludges before their disposal. Microbial masses 
are also highly non-Newtonian fluids and in several processes in 
bio- techno logy these are prepared in fixed bed reactors. Besides, 
many industrial processes involving non-Newtonian fluids are 
carried out in packed and fluidized beds to improve the heat and 
mass transfer rates. In view of these overwhelming practical 
applications cited above numerous studies have been reported on 
this subject. Reviews of published literature in the area of 
non-Newtonian flow through packed beds and porous media have been 
presented by Savins (1969), Kemblowski and Michniewicz (1979), 
Kumar et al . (1981), etc. It is appropriate to mention here that 
the motion of particles in non-Newtonian fluid is sign! f icantly 
different from its Newtonian counterpart, since non-Newtonian 
fluids exhibit a wide spectrum of rheological complexities 
Including shear dependent viscosity, visco-elasticity, time 
dependency, coupling of rheological behaviour with flow geometry 
etc. Therefore, analysing the non-Newtonian flow la a much more 
complex and challenging task than that for the Newtonian system 
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1.2 GENERAL TRENDS OF RESEARCH ON FLOW PAST A MULTI -PARTI CLE 

ASSEMBLAGEi 

Despite the wide spread occurrence and pragmatic importance 
of multi-particle fluid flow systems, the progress made towards 
the modelling of transport phenomena associated with these systems 
has been rather slow and Inadequate, especially for non-Newtonian 
media. The available literature on this subject is inundated with 
a great number of experimental and semi-theoretical correlations. 
Beginning with the celebrated work of Darcy (1856), numerous 
experimental and theoretical studies have been attempted to 
determine the characteristics of fluid flow through packed beds in 
general and pressure drop across the bed In particular. 
Similarly, the chief objective of research studies on 
sedimentation and fluidization processes has been the evaluation 
of the settling velocity/ fluidization velocity for various 
concentrations of solid particles. However, in all these studies, 
the detailed flow field around a typical particle in the 
assemblage is very scarcely reported beyond the creeping flow 
regime, even for Newtonian fluid flow. For non-Newtonian media , 
detailed flow field results have not been obtained so far even for 
creeping flow regime. There is a strong need to fill In the 
existing lacuna on the details of fluid particle interaction for 
both Newtonian and non-Newtonian media. 

In order to model the multi-particle effects on flow through 
a bed of particles, two approaches have been commonly adopted, 
namely (a) Capillary tube bundle model and (b) unit cell model. A 
few other approaches such as the continuum model and the submerged 
object model have also been In use albeit to a lesser extent. The 
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main obstacle to be overcome while modelling such flow problems is 
that of providing a satisfactory treatment of the 
particle-particle interactions Each of the available models in 
the literature accounts for the multi-particle effects in a 
different manner and there are advantages and disadvantages 
associated with each approach . 

T n "Capillary tube" bundle model , the interstitial space 
between the particles is visualised as a bundle of capillaries , 
whose lengths and diameters are such that it presents the same 
resistance to fluid flow as the actual interstices in the real 
multi-particle system. The starting point for the modelling is 
the well known Hagen- Poiseuille equation ( or its equivalent for 
Non-Newtonian flow) in conjunction with the concept of hydraulic 
radius . The success of this model with Newtonian fluids has 1 ed 
to its ext ension to non- Newtonian flow by using a suitable 
rheological equation to represent fluid behaviour . However, the 
success of this approach in both Newtonian and non-Newtonian 
systems has been somewhat surprising because it presents an over 
simplified picture of the real flow configuration and the concept 
of "mean hydraulic radius ” has little theoretical justification. 

In cell model , each part icl e , (assumed to be uni f ormly- 
spaced in the assemblage) is enveloped by a hypothetical fluid 
cell which represents the inter-particle interactions . The radius 
of the envelope is such that the voidage of each cell is equal to 
the statistically-averaged overall voidage of the assemblage. In 
the free surface cell model proposed by Happel (1958), the shear 
stress and the mass flux vanish on the cell surface. In a similar 
model independently proposed by Kuwabara (1959), zero vorticity is 
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prescribed on the cell surface Several investigators have 
pointed out the conceptual dl ££ lculti e. arising In the 
zero-vorticity model due to prescription of vorticity at the cell 
surface. The free surface cell model has been successfully 
applied to predict the pressure drop incurred during the flow of 
Newtonian and purely viscous type non-Newtonian fluids through 
fixed beds. It has also been employed with success in the 
analysis of the relative motion between a swarm of gas bubbles and 
a Newtonian or non-Newtonian continuous phase. However, all the 
aforementioned applications are limited to the creeping flow 
regime only. Very little theoretical work has been reported 
regarding the suitability of free surface cell model outside the 
creeping flow regime. 

1.3 SCOPE OF PRESENT WORK 

It appears from a close scrutiny of the existing literature 
that not only very little is known about the flow of shear 
thinning fluids through particle assemblages, but most of the 
earlier investigators have concentrated mainly on the overall 
macroscopic behaviour. No insight has been provided as to the 
causes of the observed macroscopic trends and no efforts have been 
made to relate these trends with the struc.ure of the flow field 
prevailing around a typical particle In an assemblage In fact, 
beyond the creeping flow regime even macroscopic quantities such 
as drag coefficients, rate of hindered settling etc. are not 
known. Bearing in mind that momentum transport serves as a 
precursor to the modelling of heat and mass transport, it is felt 
that there is a need to launch a systematic study of the detailed 
flow field structure (Pressure, velocity, streamlines and stress 
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profiles) over a wide range of physical, kinematic and rheological 
conditions The data accrued from such a study would facilitate 
the meaningful interpretation of experimental results and 
empirical correlations concerning heat, mass and momentum 
transport. For instance, the answer to the question such as when 
and why viscous drag becomes dominant over the pressure drag can 
come only from the knowledge of detailed flow field around a 
particle The influence of various factors such as voidage, 
non-Newtonian behaviour index and Reynolds number on important 
macroscopic quantities such as drag coefficients, fluidization 
velocity or settling velocity of suspension can be elucidated. 
Moreover, the flow configuration of a multi-particle assemblage 
provides an opportunity to examine the applicability of various 
constitutive relations containing parameters determined entirely 
from viscometric flow. However, it deserves to be mentioned at 
this stage that the results accomplished through such a 
theoretical study are likely to show some deviations from reality 
for the following reasons. 

1) Idealisations and approximations in modelling the 

hydrodynamic interactions between particles. 

2) Uncertainty concerning the applicability of the rheological 

constitutive equation which has been developed from 

viscometric flow data. 

3) Variance in the experimental data of different Investigators 
arising primarily from the improper characterization of 
non-Newtonian fluid behaviour. 

Uith regards to Newtonian fluid flow also, it is evident that 
very little is known about the usefulness and applicability of the 
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cell model approach In general and of free surface cell model in 
particular, to the intermediate Reynolds number regime. In view 
of the wider acceptance of Happel's cell model in a large variety 
of problems, an improved understanding in this direction is 
desired . 

In the present study, comprehensive numerical results for 
steady flow of both Newtonian and non-Newtonian (Power law type) 
fluid flow relative to an assemblage of spherical particles have 
been obtained, by the use of the free surface cell model. The 
Intermediate Reynolds number regime has been analysed for the 
Newtonian fluid, while the creeping as well as intermediate 
Reynolds number flow regimes have been considered for the power 
law type shear thinning fluid. Detailed flow field results have 
been presented and the influence of the flow field structure 
around the particle on macroscopic quantities of Interest such as 
drag force have been elucidated. A wide range of physical, 
kinematic and rheological conditions have been considered and the 
effects of problem parameters upon both microscopic and 
macroscopic flow quantities are discussed in detail. 



CHAPTER I! 


LITERATURE SURVEY 

A brief literature survey comprising a general overview of 
flow past multi-particle assemblages and a detailed discussion on 
modelling methodologies in particular, is presented here. Keeping 
in view the vastness of the subject, only pertinent literature has 
been referred to and discussed. 

2.1 OVERVIEW OF GENERAL RESEARCH ACTIVITY ON FLOW THROUGH MULTI - 

PARTICLE SYSTEMS! 

Ever since the pioneering work of Darcy (1856) the 
characteristics of fluid flow through packed beds and the 
determination of pressure drop across the bed have been the 
subject of numerous experimental and theoretical studies. The 
large number of variables involved makes exact treatment of the 
problem exceedingly difficult. In most of the available pressure 
drop correlations, the following factors have been considered: 
fractional void volume, type of fluid flow, particle shape, 
roughness, size distribution and manner of packing, etc. The 

effects of side walls and the bed support have also been reported. 
However, the theoretical treatments are usually confined to the 
simplest system involving uniformly sized spherical particles 
without wall and end effects. The findings of such studies 
pertaining to Newtonian fluid flow have been summarized in several 
excellent monographs on the topic (Carman, 1956; Scheidegger, 
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1960; Collins, 1961 ; Bear, 1972 , Dullien, 1979 ; Greenkorn, 1983; 
etc.) Reviews of literature that has evolved in the area of 
non-Newtonian flow through packed beds and porous media have been 
presented by Savins (1969), Kemblowski and Michniewj.cz (1979), 
Kumar et al . (1981), etc 

Research on the sedimentation process has been mainly 
concerned with the characteristics of the particle settling 
velocity. Theoretical and experimental analyses relating the 
particle settling velocity to the concentration of solids have 
been developed in the dilute limit (Burgers, 1941, 1942; 
Batchelor, 1972; Davis and Birdsel 1 , 1988) The intermediate and 
high concentration regimes are more important in practice but less 
amenable to rigorous theoretical analysis. They have been 
investigated experimentally by Richardson and Zaki, 1954a; Gar side 
and Al-Dibouni , 1977; Mor itomi et al . , 1986; etc. These authors 
have developed reliable empirical correlations for predicting the 
settling velocity in a container with closed bottom, as a function 
of solid concentration. Amongst the numerous equations available, 
the correlation of Richardson and Zaki (1954) has been most widely 
adopted due to its simplicity. A comprehensive review of recent 
advances in the area of sedimentation was presented by Davis and 
Acrivos (1985). A comparative study of several models for 
gravity-driven sedimentation of bi-disperse and Tri-disperse 
concentrated suspensions was recently reported by Al-Naffa and 
Selim (1989). 

Similar to the sedimentation process, the numerous practical 
applications of fluidization have stimulated theoretical and 





experimental investigations ( Joshi , 1983; Richardson, 1971) on 

this subject. Correlations giving the relationship between the 
sedimentation velocity and particulate concentration (Richardson 
and Zaki , 1954 ; Garside and Al-Dibouni , 1977 , Moritomi et al . , 

1986) have been found to be appl i cabl e to the fluidization process 
also, with the replacement of settling velocity by fluidization 
velocity. Most of the previous work concerning the modelling of 
these aspects of fluidized beds has been summarized by several 
authors (Leva , 1957 ; Kunni and Levenspi el , 19 6 9, Couderc , 1986) . 

Research on the hydrodynamics of non-Newtonian fluid particle 
syst ems has been focussed hitherto chiefly on packed beds . Even , 
the two most important design parameters for a fluidized bed, 
namely , fluidization velocity and bed expansion characteristics 
(Yu et al . , 1968 ; Mishra et al , 1975; Kumar and Upadhayay , 1981 ; 

Tonini et al . , 1981 ; Machac et al . , 1986 ) have received very 

little attention and more importantly the emerging picture is 
rather incoherent . Similarly, although many suspensions of 
particles have been known to exhibit non-Newtonian behaviour, 
investigations of sedimentation in non-Newtonian fluid media 
(Balkrishna et al . , 1971 ; Allen and Uhlherr , 1989 ) are limited 

The only theoretical study concerning the influence of 
non-Newtonian behaviour on sedimentation velocity of particles is 
that of Kawase and Ulbrecht (1981c). It was found in this study 
that pseudo-plasticity decreases the sedimentation velocity and 
its reduction is pronounced at large voidage. A vast gap still 
continues to exist in the understanding of the hydrodynamics of 


non-Newtonian flow through a multi-particle assemblage. 
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2.2 MODELLING TECHNIQUES 

Modelling the fluid flow past an assemblage of particles has 
remained a formidable task, because of the large number of factors 
Involved such as the nature of fluid, fractional void volume, 
hydrodynamic interaction between neighbouring particles, wall and 
end effects, etc The difficulty encountered in isolating the 
influence of each variable has marred the achievement of a clear 
understanding of the problem. However, Innumerable efforts have 
gone into analysing the various aspects of multiparticle fluid 
flow systems 

Some of the popular approaches employed in the modelling of 
viscous fluid flow past an assemblage of particles are the Darcy's 
law approach. The Geometric Pore model approach, the Submerged 
Object approach, Continuum approach and the Cell model approach. 
The last two approaches, namely, the continuum and the cell 
models, are completely theoretical in the sense that they employ 
only the basic principles of fluid dynamics to account for the 
hydrodynamic interaction between particles. Rest of the models, 
on the other hand, are based fully or partially upon empirical 
reasoning. 

Considerable efforts expended in analysing the non-Newtonian 
fluid flow through multi-particle systems involve an extension of 
the modelling techniques which were found successful for Newtonian 
fluids, with the incorporation of an appropriate constitutive 
equation to describe the non-Newtonian fluid behaviour. Detailed 
survey of the relevant literature shows that the effect of shear 
thinning (or pseudo-plasticity) on pressure drop and drag 
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coefficient has received considerable attention In these studies. 
There are only a few cases (Marshall and Metzner , 1967; Park, 
1972, Machac and Dolejes, 1982) where the effects of 
viscoelasticity have been accounted for, during the interpretation 
and correlation of results. However, a completely satisfactory 
picture has not yet emerged. 

Ulth regard to the modelling technique for non-Newtonian flow 
through particle-assemblages, three approaches have been generally 
adopted, namely, the geometric pore model, the Darcy’s law 
approach and the cell model approach. Among other available 
approaches, the elegant dimensional analysis of Slattery (1967) 
for visco-elastic fluids and the correlation procedure by Hassell 
and Bondi (1965) are noteworthy. Savins (1969) in his excellent 
review paper on this subject, has also suggested a method using 
the capillary model-hydraulic radius concept which does not 
involve any rheological equation. 

2 . 2 . 1 Geometric Pore Model 

A geometric pore model which has been extensively used to 
idealize the flow through the packed bed is the "Capillary tube 
bundle” approximation In capillary model, the space between 
particles is represented in terms of a bundle of tortuous 
capillaries, the length and diameter of which is determined such 
that the same resistance to flow is offered as the actual 
interstices of the multi-particle system. The starting point for 
this model is the Hagen-Poiseuil le equation. The moat well- 
known representations of this class of models are due to Blake, 
Blake-Kozeny , Kozeny-Carman and Ergun. In Blake’s model, the bed 
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le envisaged as a bundle of straight tubes of complicated 
cross-section The Blake-Kozeny model portrays the bed or porous 
medium as a bundle of tangled tubes with an effective length 
greater than the actual length of the bed by a tortuosity factor 
In the Kozeny- Carman model, the effect of tortuosity on the mean 
velocity of flow is also considered by taking the mean velocity as 
the product of Dupuit ( 1863 ) velocity and the tortuosity factor . 
According to this model , the final expression for the modified 
friction factor , f^, in terms of the modi f i ed Reynolds number, 

Re^ takes the simple form 

£ bk - *57— C2 15 

bk 


The Kozeny-Car man model may be considered to provide a more 
accurate representation of flow through the bed , comprising of 
tortuous tubes . 

All the above discussed capi 1 lary models are applicable for 
laminar flow regime only. The correlation for pressure drop 
through a fixed bed in turbulent flow regime was first developed 
by Burke and Plummer (1928). Subsequently , Ergun (1952) combined 
this correlation with that of Kozeny- Carman to provide a single 
expression valid for the entire flow range . The Ergun equation 
which essentially considers the total pressure drop as the sum of 
viscous and inertial pressure losses through the bed is expressed 
in the form. 


f 


bk 



+ 1.75 


( 2 . 2 ) 


The success of capi 1 lary models for Newtonian fluid flows has 


led to their extension to non-Newtonian flows also by many 
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investigators The earliest attempts in this direction were made 
by Christopher and Middleman (1965), who used the power law fluid 
model; and Sadowski and Bird (1965), who in a parallel treatment, 
used the Ellis fluid model to study the creeping flow. These 
authors related the flow rate and pressure drop by suitably 
modifying the Blake-Kozeny equation. Since then , numerous 
investigators have followed suit and studied a wide variety of 
fluid models (Gaitonde and Middleman , 1967; Kemblowski and Mertl , 
1974 , Brea et al . 1976 ; Kumar and Upadhyay , 1981; Kozicki and Tiu, 
1988). Most of these studies are restricted to the creeping flow 
regime, save those of Brea et al . (1976) and Hanna et al . (1977) 
who obtained equations covering the laminar as well as turbulent 
flow regimes . Srinivas and Chhabra (1990) recently made 
exhaustive comparisons of the predictability of methods available 
in literature f or the prediction of pressure loss during the flow 
of non-Newtonian fluids through packed beds . In this study, they 
used their own experimental data and those of others ' available in 
literature (Sadowski and Bird , 1965; Park et al . , 1975 ) . The 
range of power law index covered in this study was 0.38 to 1 . 0 and 
the fixed bed porosity varied from 0.369 to 0.447. It was 
concluded that almost all the correlations overpredict the value 
of pressure drop for a given flow condition. The correlation 
proposed by Dharmadhikari and Kale (1985) was found to provide the 
closest agreement amongst all the results examined. 

Despite the popular use of capillary models and numerous 
experimental results available in literature, there is still no 
agreement on the value of the constant A occurring in the 
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Kozeny- Carman equation. For Newtonian fluids, the values 
suggested range from 115 to 200, but the most commonly used are 
the values of 180 and 150 as determined by Carman (1937) and Ergun 


(1952) respectively 

Kemblowski and 

Michniewicz 

(1979) 

have 

r e-analyzed 

most of 

the experimental 

data 

in 

literature 

and 

suggested a 

value of 

180 for Newtonian 

and 

150 

for power 

law 


fluids. However, based on their own analysis, they have suggested 
a value of 180 for power law fluids, which was subsequently 
confirmed also for the flow of Carreau fluids (Kemblowski and 
Michniewicz, 1981). Apart from the afore-mentioned inconsistency, 
the capillary models have been criticized by several investigators 
(Barnea and Mednick, 1978; Molerus, 1980) for the following 
reasons • 

(a) The tube bundle model represents an over-simplification of 

the real porous medium and the concept of mean hydraulic 
radius has little theoretical justification for laminar flow 
in a non-circular channel. Sheffield and Metzner (1976) and 
many others found the capillary model to be conceptually weak 
in that it does not take into account the converging- 
diverging nature of flow path and the inter-connection 
between pores. In case of non-Newtonian fluids this 
deficiency appears to be more serious because the 

viscoelastic effects are manifested in a converging-diverging 
flow path. 

(b) In laminar tube f 1 ow the inertial ef f ects are obs erved only 
with the onset of turbulence beyond a critical Reynolds 
number. In the case of flow through particle assemblages, 
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due to frequent deflection of fluid flow, inertial effects 
begin to dominate long before the onset of turbulence. In 

this regard, the pipe flow similarity with the flow in a 
packed bed does not exist 

(c) Durst et al (1987) have pointed out that the pressure drop 
of porous media flow is only to a small extent due to the 
shear force term which is usually employed in the capillary 
model. For more correct derivations, additional terms have 
to be taken into account since the flow field is subjected to 
both shear and elongational strain. 

In order to partly obviate some of these deficiencies, 
complex pore models have been proposed (Batra et al . , 1970 
Dullien, 1975; Scheidegger, 1974; Payatakes and Neira, 1977; 
Sheffield and Hetzner , 1976; Deiber and Schowalter, 1979; 
Kanel lopoulos , 1985). Although such more sophisticated pore 

models may be useful in accounting for some of the second order 
effects, they do not offer any measurable improvement (Cohen and 
Hetzner , 1981) concerning the prediction of r elt ionships between 
the pressure drop and the flow rate. Therefore, classical 
capillary models continue to remain popular for describing flow 
through packed beds. 

3m 2. 2 Darcy* s Modal 

Most of the available analytical work on transport processes 
through porous media are based on Darcy's law. This law was 
established in 1856 by Darcy empirically, and since then has been 
verified experimentally by numerous investigators for Newtonian 
fluid motion in porous media at low flow rates. It states that 
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for an isothermal fluid moving with a slow steady velocity V 
under a pressure gradient , through a homogeneous and isotropic 

porous bed 

V = - £ (2 3) 

Yi ax 

where, y> is the fluid viscosity and C is an empirical constant 
called the permeability of the medium. Originally, the Darcy’s 
law was established to obtain the average or apparent flow rates 
of fluids through porous media in one-dimensional flows Later 

researchers, assuming the fluid velocity to be a continuous 
function of space coordinates, have given Darcy’s law for two or 
three dimensional flow fields. The extension of Darcy’s law to 
multi-dimensional flow including the body force term is expressed 
as , 

V = £ [ - Vp + f ] (2.4) 

y? g 

where, f is the body force vector given by the weight of the 
& 

fluid per unit volume. It is to be noted, however, that the 

Darcy’s law is established empirically and it can not be derived 
analytically by performing the momentum balance on a small element 
of porous medium. It is merely a constitutive relation which does 
not yield much information about the permeability of the porous 
medium. The drawbacks of Darcy's law are that it does not satisfy 
the no-slip boundary conditions on impermeable side walls 
confining the porous medium and also It becomes inapplicable if 
the permeability of the porous medium is high. In these cases, 

due consideration needs to be given to factors such as inertial 

effects in the fluid. Later generalizations of the Darcy’s law 
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have focussed attention on some of these factors. 

The inclusion of viscous shear stress terms in the flow 
equations and the consideration of a velocity field which 
satisfies no-slip at impermeable walls , have been proposed by 
Brinkman (1947). This is the so called Brinkman’s extension, 
which accounts for the flow resistance due to the fluid viscosity 
itself. Yamamoto and Iwamura (1976) proposed a flow model in 
which the effects of inertia are taken into account via convection 
terms in the momentum equations Obviously, this is an important 
correction which has to be considered at higher flow rates through 
the porous medium Yet another approach to handle the Inertial 
effects was proposed by Forchheimer (1901), which approximately 
accounts for the pressure loss due to flow separation in the 
pores, by the inclusion of an additional term in the flow 
equations. Ulth all these extensions, the flow models for porous 
media situations have more or less the same general structure as 
that of the Navi er -Stokes equations for free fluid flows. It is 
to be mentioned that despite these developments most of the 
studies available at present on porous medium applications, 
especially for non-Newtonian fluid flow, do not consider a 
generalized Darcy’s law. 

The Darcy’s model has also been extended to non-Newtonian 
flows with the help of suitable rheological equations. This 
approach has been adopted by Kozicki et al . (1967, 1973), Gregory 
and Griskey (1967), Al-Farris (1989), and some others for a 
variety of fluid models. 
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2. 2. 3 Submerged Object. Approach 

The capillary and Darcy models are reported to fail at higher 
voidage and the experimental facta are better described by the 
submerged object models in this range These models are also 
known as discrete particle models, as they involve the 
determination of pressure drop across the bed from the expression 
for drag on a single particle immersed in the fluid, by assuming 
simple drag force additivity. Hydrodynamic interactions of 
neighbouring particles which affect the drag force on each 
individual particle, are not taken into account. The only 
indirect effect of neighbouring particles that is taken into 
account is the decrease in the average cross-sectional area 
available for fluid flow. 

At a very early stage, Burke and Plummer (1928) had proposed 
a discrete particle model which regards each particle in the 
assemblage as being enveloped by a separate boundary layer. This 
approach was later expanded by Ranz (1952). Subsequently many 
other variations (Gauvin and Katta, 1973; Foscolo et al . , 1983, 
Barnea and Mizrahi, 1973; Barnea and Mednick, 1978; Zimmels, 1988) 
of the submerged object model have been reported. The method of 
determining the drag force on a typical particle in the assemblage 
differs from one model to another; however, the common feature of 
all the submerged object models is the concept of simple drag 
force additivity. 

The submerged object models are comparatively less popular 
than the capillary models and cell models, since they do not 
incorporate the particle to particle interaction 


directly. 



21 


However, analyses of non-Newtonian multiparticle systems based on 
the model are not available in the literature. 

2.2.4 Continuum Approach 

In the continuum model, the hydrodynamic forces and torques 
exerted on the particle by the fluid are computed by taking into 
account all part 1 cl e- part i cl e and particle-fluid interactions and 
then solving this complicated problem in the bounded region. One 
of the earliest exact solutions accounting for particle 
interactions in this way is that of Stimson and Jeffery (1926). 
They solved the problem of slow motion of two equal-sized spheres 
parallel to their line of centres. For the same problem Cooley 
and O’Neill (1969) gave a more comprehensive set of results for 
the cases in which spheres may be of unequal sizes. An exact 
solution for the two spheres problem has been possible primarily 
because of the use of separation of variables technique m 
bispherical coordinates. For a larger collection of spheres, or 
for a pair of non-spher ical particles, it is not possible to find 
a suitable coordinate system for this purpose. It is then 
necessary to seek an approximate solution to the problem using a 
systematic scheme of successive iterations, thereby the boundary 
value problem may be solved in principle to a degree of accuracy 
by considering boundary conditions associated with one particle at 
a time. Such schemes are provided by the "method of reflections” 
(Smoluchowskl , 1972) or the pairwise additive interaction model 
(Kim and Russel, 1985). 

The major difficulty in modelling the multi-particle 
interaction by this model arises in summing the contribution of 



different particles on the test particle. The analysis is 
complicated by the fact that in a system containing a large number 
of particles, the resulting sums or integrals are often divergent. 
One of the earliest attempts to overcome this difficulty was made 
by Hasimoto (1959), who obtained spatially periodic (Fourier 
series) fundamental solutions to the Stokes equation of motion for 
viscous flow past a periodic array of spheres This approach was 
later employed by Zick and Homey (1982) and modified by Sangani 
and Acrivas (1982) to Investigate slow flow past an array of 
spheres. Amongst the several other methods available (Hinch, 
1977) to handle the divergent nature of sums of all contributions 
on the test particles, the renormalization technique (Batchelor, 
1972; O'Brien, 1979; Feuillebois, 1984; Kim and Russel, 1985) has 
been widely used. 

In summary, though it has been possible in principle, to 
formulate continuum models for most multi-particle flow problems, 
the solutions are generally too complicated to be useful or 
sometimes are even impossible to obtain. In general, the 
applicability of this model has so far been limited to dilute 
multiparticle systems and low Reynolds number flows of Newtonian 
fluids. Investigations of non-Newtonian fluid flows past an 
assemblage of particles using continuum approach is hardly 
available in literature. 

a. 2. 5 Cell Model 

Amongst the various models to describe the hydrodynamic 
interactions between the particles of a multi-particle system, the 
cell model has received extensive attention, probably because of 
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the relative ease with which significant results can be obtained. 
In its formulation, the difficult many body problem is replaced by 
a simple and conceptually more amenable situation Involving one 
particle. Wall, entry and exit effects are usually neglected. It 
is then of great application in concentrated assemblages, where 
the effects of container walls will not be important. The 
assembly of particles in the fluid is usually assumed to be 
uniform and each sphere is fixed in space The interaction of a 
particular sphere with its neighbours is modelled by considering a 
hypothetical fluid envelope around that sphere such that the 
fractional void volume in the cell is identical to that of the 
entire assemblage. The two most extensively used cell models are 
the Happel’s free surface model (1958) and the Kuwabara * s zero 
vorticity cell model (1959). However, many other cell models have 
also been reported in the literature. 

A somewhat similar model for study of the velocity of steady 
fall of spherical particles through a fluid medium was employed by 
Cunnigham (1910), who assumed that each particle In a cloud would 
be effectively limited to motion within a concentric mass of 
fluid. However, he assumed that the boundary of the outside fluid 
envelope was solid; this corresponded in someway to the surfaces 
of the surrounding spheres present In the cloud. As commented by 
Happel (1958), in this model the size of the spherical envelope 
must be fixed by some additional considerations. 

Uchida (1949) proposed a model for viscous flow through 
packed beds. He assumed that the bed is divided Into numerous 
centres and the particle is placed in the centre of each cube. 
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The boundary value problem of satisfying the appropriate boundary 
conditions simultaneously on spherical and cubical surfaces proved 
Intractable and consequently, Uchida simplified it. The 
prediction of Uchida for the pressure drop through the cubic 
assemblage of particles does not agree well with experimental 
results of Happel and Epstein (1954) and the fixed bed data 
represented by Kozeny-Carman equation (Carman, 1937). 

Richardson and Zaki (1954b) developed a semi-theoretical cell 
model for sedimenting spherical particles. The particles in the 
suspension were arrayed in a hexagonal pattern in horizontal 
planes. Two models were proposed In the first, the vertical 
layers had the same spacing as that between the particles in the 
horizontal plane and in the second, the vertical spacing was such 
that adjacent horizontal layers of spheres were in contact. In 
order to simplify the boundary conditions, it was assumed that 
each particle was surrounded by a cylinder of fluid with the same 
cross sectional area as the hexagonal prism. The appropriate 
boundary conditions for symmetric flow with respect to adjacent 
particles and zero vorticity at the particle surface, were 
prescribed for the system and solutions were obtained. In the 
development, the theoretical treatment was not rigorous and many 


simplifying assumptions 

were 

made . 

The 

predictions 

of 

sedimentation 

velocity 

from 

this 

analysis 

for the 

two 

conf igurat ions 

were found 

not 

to be 

applicable to a dilute 


concentration of particles, although configuration II gives better 
agreement of the two with the experimental data of Richardson and 
Zaki (1954a) in the intermediate particle concentration range. 
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The deviations may be attributed partly to the over-simplification 
of the pressure gradient and partly to the questionable 
assumptions (Jean and Fan, 1989) regarding the angles between the 
stream lines and the vertical axis . 

Kavaguti (1958) proposed a "cylindrical cell model" to 
estimate the sedimentation velocity of uniform spheres in the 
Stokes flow regime. The drag force on a particle exerted by the 
fluid in a cylindrical cell was evaluated, based on viscous flow 
of a fluid past a sphere in a frictionless circular pipe of 
infinite length Since Faxen’s procedure (Faxen, 1923) was used 
to approximate the velocity profile for such a flow, the result is 
valid only for a relatively dilute particle concentration. 
Furthermore, the relationship between the particle concentration 
and the ratio of the particle diameter to the radius of the 
Imaginary pipe, is based on the effective projected area of the 
pipe and particle rather than the volume, an assumption whose 
validity is rather questionable. 

Happel (1958) proposed a "spherical" cell model similar to 
that of Cunnigham (1910), and attempted to predict the bed 
expansion (or sedimentation velocity in a sedimentation process) 
and the pressure drops at low Reynolds numbers for fluidized and 
fixed beds respectively. The only difference in Happel’s cell 
model compared to that of Cunnigham (1910) is that zero shear 
stress is prescribed at the outer spherical fluid envelope. This 
assumption is probably more reasonable than that of zero fluid 
velocity on the outer sphere adopted by Cunnigham. The 
specification of zero shear stress boundary condition is motivated 
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by the reasoning that the cells are of non-interacting type, 
thereby requiring no energy exchange with the surroundings. The 
entire disturbance due to each particle is thus confined to the 
cell of the fluid with which it is associated Within the bounded 
region of the cell, the Navier-Stokes equations after omitting the 
inertia terms, describe the fluid motion. The Happel cell model, 
also known as the free surface cell model , has been successful in 
predicting the pressure drop and hindered settling characteristics 
of concentrated (0.3 < « i 0.6) systems of spheres. It reduces to 
Stokes law as a limiting case at infinite dilution, but fails to 
predict the hindered settling behaviour of dilute suspensions. 
For the fluidization or sedimentation cases, it shows poor 
agreement with the experimental data of Steinour (1944), Virschoor 
( 1951 ), Richardson and Zaki (1954a) and Hanratty and Bandukwala 
( 1957 ), but shows reasonable agreement with the data of Adlei 
( 1958 ). It is important to recognise that the type of flow 
distributor used in the experiment exerts a significant influence 
on the flow field (Jean and Fan , 1989). 

Theoretical predictions based on the Happel ’ s model agree 
well with the experimental data for heat and mass transfer in 
multiparticle systems in the creeping flow region (Pfeffer, 1964; 
Pfeffer and Happel, 1964). The Happel ’ s analysis for solid 
particles has been extended to analyse swarms of drops and bubbles 
in the presence or absence of surfactant impurities (Gal-or and 
Uaslo , 1968). Subsequently, Yaron and Gal-or (1971, 1972) used 
the hydrodynamic solution of Gal-or and Uaslo (1968) to 
investigate heat and mass transfer processes from side-distributed 
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drops, bubbles or solid particles and later on to calculate the 
effective viscosity of concentrated suspensions and emulsions. 

Another spherical cell model was proposed by Kuvabara (1959), 
who used a zero vorticity assumption instead of the zero shear 
stress assumption on the outer spherical envelope in otherwise 
identical flow conditions. He calculated the drag in the creeping 
flow regime by integrating the dissipation energy function. 
Several authors (Happel and Brenner, 1965; El-Kaissy and Homsy, 
1973) have criticized the zero vorticity spherical cell model for 
the inconsistency arising due to the fact that shear stress does 
not vanish at the outer surface; this is tantamount to the 
exchange of energy between neighboring cells, which violates the 
basic promise regarding the non-interacting nature of cells. 
Also, the zero vorticity cell model results in larger energy 
dissipation within the fluid envelope as compared to the free 
surface cell model of Happel (1958); the increased dissipation 
presumably is due to the additional work done by the stresses 
existing at the outer fluid envelope. As a consequence, 
Kuwabara ' s solution yields a stronger concentration dependence of 
the drag on a typical particle, experienced by uniformlly 
distributed spheres in viscous flow at low Reynolds number. 

Most of the applications of cell models in describing the 
transport processes in multi-particle assemblages have been 
limited to creeping flow conditions . Very little work has been 
reported outside the creeping flow regime. Leclair and Hamielec 
(1968) solved the complete Navier-Stokes equations for the zero 
vorticity cell model using the finite difference method and 
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reported the values of drag coefficient as a function of bed 
voidage and Reynolds number, encompassing wide range of conditions 
(0.408 S « * 0 9; 1 £ Re £ 500 ). It is not at all clear whether 
the choice of zero vorticity cell model was based on its merit or 
simply because of its ease in Implementation for vort icity-str earn 
function approach. In an attempt to explore the extent to which 
the spherical cell models may be extended to higher Reynolds 
numbers, El-Kaissy and Homsy ( 1973 ) carried out a regular 
perturbation analysis for both free surface and zero vorticity 
cell models and have provided analytical expressions for the 
friction factor as a function of bed voidage (0.3 S « ^ 0.6) and 
Reynolds number. However, the perturbation scheme has limited 
applicability beyond the creeping flow regime in the sense that as 
the interparticle distance increases (i.e. the voidage increases), 
the maximum value of Reynolds number upto which their analysis is 
valid decreases Based on the comparison between theory and 
available experimental data, El-Kaissy and Homsy found that the 
free surface cell model yields consistently better results than 
the zero vorticity cell model. However, for Re > 20, both models 
show appreciable departures from experimental measurements. 

Nlshimura and Ishii (1980a, b) employed the free surface cell 
model together with the Karman-Pohlhausen integral method to 
examine the momentum and mass transfer boundary layers in particle 
assemblages at high Reynolds numbers. Using fourth order 
polynomial expansions these authors reported extensive results on 
the variation of friction factor and Sherwood number as functions 
of bed voidage and Reynolds number and claimed their results to be 
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applicable in the range 10 £ Re £ 10 . Intuitively, it ia 
unlikely that the boundary layer approximation will be valid at 
values of Reynolds number as low as 10. Another somewhat unusual 
feature of their study is that the form drag is shown to be 
independent of the value of Reynolds number ; also their results in 
the limit of voidage approaching unity differs from the value of 
drag coefficient available In literature by a factor of 2 or more. 
In a subsequent paper Fukuchi and Ishii (1982) employed a 5th 
degree polynomial and re-investigated the same problem However, 
qualitatively similar results have been reported in this study 
also 

In addition to the aforementioned widely used cell models, 
many other cell models have been reported but their applicability 
has been examined only for a narrow range of applications. Gunn 
and Malik (1966) proposed a model for fluidized beds and 
sedimenting systems based upon arrays of particles where periodic 
conditions are established throughout the flow field. This model 
postulates the enclosing boundary to be frictionless; therefore, 
the total pressure drop across the cell can be equated to particle 
drag. As the mathematical difficulty of solving the equations 
subject to the boundary condition was formidable, Gunn and Malik 
(1966) presented only some general properties of solution in 
conjunction with the experimental results. 

Like the sphere-in-sphere configuration of spherical cell 
models, sphere- in-cy Under configuration leading to what is 
commonly called as "cylindrical cell model" have been proposed by 
several investigators. Such models were earlier proposed by 
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Happel and Ast (I960), Kavaguti (1958) and Uakiya (1953). The 
Happel and Ast model was developed to represent the process of 
sedimentation The key concept in this model is the postulate of 
zero shear at the cyl inder walls Within the bounded region, the 
Stokes equation for viscous flow was solved by using the method of 
Haberman to obtain the drag on the particle. It is important to 
add here that all these cell models apply in creeping flow regime 
only Recently, a cylindrical cell model for the hydrodynamics of 
particle assemblages at intermediate Reynolds number was reported 
by Tal and Sirignano (1982). These authors suggested that their 
model can be applied to study transport process in fluidized and 
packed beds as well as In predicting the rheological properties of 
suspensions . 

The success of Happel ' s free surface cell models (1958) In 
the Newtonian case has led to its natural extension to 
non-Newtonian fluids. For instance, Mohan and Raghuraman (1976a, 
1976b) extended its applicability to power law and Ellis model 
type fluids to calculate the drag coefficient Later on, an 
analogous treatment for the Carreau model fluids was presented by 
Chhabra and Raghuraman (1984). Qualitatively, for power law type 
fluids the analysis of Mohan and Raghuraman (1976a) predicts an 
Increase in the value of drag coefficient above its Newtonian 
value, whereas, if a fluid model containing a zero shear viscosity 
is used to depict the shear thinning behaviour, a reduction in the 
drag Is predicted below the Newtonian value (Mohan and Raghuraman, 
1976b; Chhabra and Raghuraman, 1984). It needs to be emphasised 
here that there is no inconsistency in this Whether the drag is 
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augmented or suppressed with reference to the corresponding 
Newtonian value is only a matter of how the Reynolds number is 
defined However, a notable feature of all these analyses (which 
are approximate only) is that the two bounds diverge with the 
increasing degree of shear thinning characteristics In the 
absence of any definitive information, the use of arithmetic 
average of two bounds has been generally suggest-ed, In another 
study, Kawase and Ulbrecht (1981a) linearised the governing 
equations using the approach earlier developed by Hirose and Moo 
Young (1969) and obtained a closed form analytical expression for 
the total drag force experienced by a particle assemblage in 
relative motion with respect to power law liquids This analysis 
also pr edicts an incr ease in the drag foice over its Newtonian 
value. Owing to the approximate nature of their treatment and 
associated assumptions , the results of Kawase and Ulbrecht 
(1981a) are likely to be applicable only for mild to moderate 
shearthinning fluids through multi-particle assemblages. The only 
other study of the flow of shear thinning fluids through 
multi-particle assemblages is that of Hua and Ishii (1981). In 
this study, high Reynolds number flow of power law liquids through 
particle assemblage was simulated by employing the assumption of 
boundary layer flow in Happel ’ s free surface cell model. 

A close scrutiny of the above cited literature reveals that 
results are limited on the detailed flow field structure around a 
typical particle in an assemblage and explanations for the 
observed drag phenomena are still wanting , except for the creeping 
flow of Newtonian fluids. The performance of the Happel cell 



32 


model in the intermediate Reynolds 
ascertained As already described 
the Intent of the present study to 
in detail 


number regime is also 
at the end of Chapter- 1 
analyse some of these 


to be 
it is 
aspects 



CHAPTER II! 


MATHEMATICAL MODELLING 


The flow of a generalized Newtonian fluid relative to an 
assemblage of spherical particles has been considered. Modelling 
of such mult ipart icl e fluid flow systems rests upon the 
application of the principles of conservation of mass and 
momentum, in conjunction with a constitutive equation for the 
rheological behaviour of the fluid and a model for the 
representation of the hydrodynamic interactions between the 
particles, wall effects and end effects. The choice of a 
Generalized Newtonian fluid facilitates the use of the same 
formulation for both Newtonian as well as inelastic viscous 
non-Newtonian fluids. However, before embarking upon the actual 
formulation of the problem, it is instructive and desirable to 
present a terse discussion of the rheological equation of state 
for a non-Newtonian fluid. 

3.1 GEOLOGICAL EQUATION OF STATE 

One of the simplest forms for the constitutive relationship 
between the stress tensor and the rate of deformation tensor at a 
point in a fluid is the Newton’s law of viscosity. For an 
incompressible fluid, the stress tensor can be written as 

a 


P6 + r?A 


(3.1) 



34 


where the rate of deformation tensor is given by 

A = VV_ + (<?VJ T (3.2) 

For a Newtonian fluid, the viscosity depends only on temperature 
and pressure, and the fluids obeying equations (3.1) and (3.2) 
with a viscosity variation described by n = r>( p, T ) are called 
Newtonian fluids All gases and low molecular weight liquids are 
known to display Newtonian behaviour. 

It is well known that structurally complex fluids such as 
polymer melts and their solutions, suspensions and emulsions, etc. 
exhibit complex rheological behaviour which can not be described 
by the Newton’s law of viscosity and hence these fluids are 
accordingly called non-Newtonian fluids. From an engineering view 
point, perhaps the most Important non-Newtonian effect is the 
shear rate dependence of viscosity. Most non-Newtonian fluids 
display shearthinning or pseudo-plastic behaviour, that is, their 
apparent viscosity decreases with Increasing shear rate. Other 
non-Newtonian effects include the existence of unequal normal 
stresses in steady shear flow and complex transient response in 
unsteady shear flow. Some examples of the former are the 
"Ueissenberg effect" (Bird, Armstrong and Hassager, 1977) in which 
a macromol ecular fluid climbs a rotating rod; the "hole pressure 
error” Inherent in pressure transducers that are not flush mounted 
(Higashltani and Pritchard, 1972; HIgashitani and Lodge, 1975) and 
the reversal of the direction of secondary flow when a disc 
rotates a liquid in a beaker. Unsteadyy state responses, on the 
other hand, are mainly observed In dynamic testing, stress 
relaxation, creep, recoil (Ferry, 1970) and in the other time 
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dependent experiments such as the start up of a shear flow, etc. 

All the examples cited above are confined to shearing flows. 
Additional non-Newtonian effects which are encountered in 
non-shearing flows are the die-swell (Lodge, 1964); development of 
a toroidal vortex in the inlet flow to a die (Tordella, 1957; 
Gleskus, 1968; Metzner , Uebler and Chan Man Fong, 1968); the 
bizarre "Kaye Phenomenon", where a falling stream of polymer 
solution impinging on the puddle of the same fluid "bounces back” 
from the liquid surface (Lodge, 1964). Bird et al . (1977) have 
presented all excellent discussion of many of these pohenomena in 
their classic treatise. 

In order to describe these non-Newtonian effects 
quantitatively Eqn . (3.1) must be replaced by an appropriate 
rheological equation of state. The earliest attempts involved 
simple adhoc empiricism to describe some very limited class of 
such flows. In this catagory the models of "Generalised Newtonian 
fluid and the "Linear viscoelastic fluid" were developed to 
describe shearthinning behaviour and the transient phenomena in 
small scale deformations, respectively Further, attempts were 
made to devise empirical "non-linear viscoelastic models" that 
could describe and interrelate, at least qualitatively, a wide 
variety of phenomena. Finally more complex Integral expansions of 
basic memory functions were evolved which provide a better 
perspective of the aforementioned effects. However, it is also 
appropriate to add here that failures encountered in many careful 
numerical studies in recent years, have raised questions about the 
utility of some of these (Rallison and Hinch , 1988) constitutive 
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equations for describing the viscoelastic effects. Since the 
present work is concerned with the hydrodynamics of generalised 
Newtonian fluids (GNF ) in multiparticle assemblages, it is 
appropriate to present a brief discussion on the rheological 
behaviour of this class of fluids. 

3.2 GENERALISED NEWTONIAN FLUIDS 

The general expression used to describe the behaviour of a 
time independent isotropic fluid is given by 


& = pS + T 


(3.3) 


t — ftA + A + T) { A . A ) 


(3.4) 


where Y) is the ” non-Newtonian viscosity”, Y) is the "cross 
viscosity” and r is the deviatori c stress tensor; ft, Y) and n c are 
functions of the three principal invariants of the rate of 
deformation tensor, A. The fluids obeying Eqn. (3.4) are known as 
Reiner-Rivlin fluids. The absence of time derivatives in the 
above equation is tantamount to excluding viscoelastic phenomena 


from the theoretical consideration. 


For incompressible fluids 


both ft and the first invariant are identically zero, and 


therefore, eqn. (3.4) reduces to 


t = Y}, A^ + tj {A .A) 


(3.5) 


where r? = r> ( 1 2 » 1 3 ) ; and 1 2 and I 3 , the second and third 

invariants of the rate of deformation tensor are given by the 
expressions 


I, = A : A 


(3.6) 


I = det A 


(3.7) 


Several investigators have questioned the usefulness of the 
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Re iner-R ivl in equation on the grounds that the available 

experimental evidence does not support the existance of 

non-vaniahing normal stress coef f i ci ent r> c without appreciable 

viscoelasticity. However, it is possible to derive the 

constitutive equations for a purely viscous (inelastic), 

incompressible fluid by considering ft and r> as zero in eqn (3.4) 

c 

and Y) being a function of only the second invariant of the rate of 
deformation tensor. Thus, 

t = >) (3.8) 

where » = r> ( 1 2 ) (3.9) 

Reiner called eqn. (3.8) as the Generalised Newtonian fluid model 
and it reduces to the Newtonian case for a constant value of r ) . 
The form of eqn. (3.9) has received further support from the fact 
that the third invariant of the rate of deformation tensor is 
identically zero for viscometric flows and it is considered to be 
unimportant in other flow situations as well (Slattery and Bird, 
1961; Tanner, 1966). Hence, eqn. (3-9) represents a good 

approximation suitable for engineering applications. 

Based on experimental measurements conducted for viscometric 
flows, the functional dependence of the appareent viscosity n on 
I ^ has been portrayed by a variety of empirical formulations. 
Consequently, there is no scarcity of such equations in the 
literature and these have been summarized by Bird (1965, 1976) and 
Skelland (1969), among others. 

3.2.1 The Newtonian Model 
For a Newtonian fluid 

r) - constant = p (3 10) 
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where p is the Newtonian viscosity . Hence , 

^ = P £. (311) 

and in steady one dimensional shear flow 

r = p e (3 12) 

where e is the shear rate 


3. 2. 2 The Power-Law Model 

According to this model, 
n-1 



the apparent viscosity is given by 

(3 13) 


where n is the power law flow behaviour index and K is the 
consistency factor, which are determined using shear stress-shear 
rate data For steady one dimensional shear flow, Eqn. (3 13) can 
be re-written as 

t = K e 11 (3.14) 


The apparant viscosity is defined from Eqn (3.14) as 

r) = K | e ( n_1 (3.15) 
Eqn. (3.15) is probably "the most criticized , most maligned and 
yet most widely used equation in all of the rheology” (Schowalter, 
1978). The main virtue of the power law model is its simplicity, 
as the extent of shear thinning behaviour is simply characterised 
by the value of n. Uhen n is unity, Eqn (3.15) reduces to the 
standard Newtonian fluid with a viscosity K. Uhen n < 1, the 
viscosity decreases with increasing shear rate (pseudoplastic or 
shearthinning behaviour); and for n > 1 the viscosity increases 
with Increasing shear rate (di latant or shear thickening 
behaviour). Most fluids tend to exhibit power law behaviour over a 
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range of shear rates. One of the severe drawbacks of the power 
lav model is that It fails to predict a constant viscosity in the 
zero shear rate region. Numerous models containing more than two 
parameters are available, which overcome some of the drawbacks of 
the power law model. These have been covered adequately in 
several references e.g see Bird (1965), Skel land (1967), Darby 
(1967), Bird et al . (1977), etc. However , it must be borne in mind 
that the introduction of each parameter in the rheological 
equation of state also makes the solutions of hydrodynamic 
problems increasingly difficult (Katsuhisa and Bit d , 1965) The 
two such models which have found widespread use in approximating 
shear stress-shear rate behaviour of a range of several substances 
are the Ellis and the Carreau models Several workers have 
employed these models for investigting the behaviour of 
non-Newtonian fluid particle systems (Mohan and Raghur aman , 1976b; 
Chhabra and Raghur aman, 1984; Kemblowski and Michniewicz , 1981; 
Gummalam , Narayan and Chhabra, 1988; Bush and Phan-thlen, 1984). 
However, there is no question that in majority of the experimental 
and theoretical investigations pertaining to viscous flow past an 
assemblage of particles, power law model has been employed. In 
the present study also, the power law type Generalised Newtonian 
fluid model has been selected for describing the fluid behaviout . 
However, the problem formulation itself has been made general so 
as to be applicable to other kinds of Generalised Newtonian fluids 


with minor modifications. 
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3.3 PARTICLE-PARTICLE INTERACTION 

Varioua approaches adopted to model the hydrodynamic 
interactions between the particles have been discussed in Chapter 
II. Though in principle, it is possible to formulate continuum 
models for multiparticle systems, the solutions of these models 
are generally too complicated to be useful In contrast to this, 
the cell models have been extensively used in modelling the 
multi-particle assemblages and indeed have yielded significant 
results. As discussed in the previous chapter, the free surface 
cell model has gained wide acceptance in predicting the pressure 
drop through packed beds (Mohan and Raghuraman , 1976a; Kawase and 
Ulbrecht, 1981; Nishimura and Ishii , 1980a, b; Chhabra and 
Raghuraman, 1984), settling characteristics of concentrated 
suspensions (Kawase and Ulbrect, 1983; Prasad et al , 1990) and 
rise and fall of drops and bubbles (Gal-Or and Uaslo , 1968; 
Bhavraju et al . , 1978; Jarzebski and Malinowski, 1986; Gumma lam et 
al , 1988) in Newtonian as well as in non-Newtonian media. Hence, 
this model is employed to describe the particle-particle 
interactions in the present work. 

The free surface cell model envisions the particle assemblage 
to consist of a number of identical unit cells, each of which 
contains a particle placed inside a fluid envelope such that the 
voidage of the cell is equal to the statistically averaged overall 
voidage of the assemblage. The cells are of non-interacting type 
and therefore, one unit cell represents the complete assemblage, 
neglecting the wall and end effects. 
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ig 3 1: Schematic representation of the flow configuration 
(Happel Cell). 



3.4. PROBLEM STATEMENT AND FORMULATION 


The steady, linear motion of a uniformly spaced assemblage of 
rigid spherical partiales through an incompressible, inelastic, 
viscous fluid has been considered. The rheological behaviour of 
the fluid is modelled using the power law equation (Eqn 3 13), 
with Newtonian fluid as a special case (n = 1) Each particle of 
radius Rp in the assemblage, is isolated from the rest of the 
particles by defining a cell boundary of radius R c concentric with 
the particle. The flow around a typical particle in its cell is 
the flow configuration analysed The schematic representation of 
the flow configuration (unit cell) along with the coordinate 
system (r, & , <&) is shown in Fig. 3 1. Axisymmetry is assumed to 
prevail within the cell about the direction of particle motion. 
The solid sphere moves in the positive X-direction with a 
velocity, As stipulated by the free surface cell model, the 
radial velocity and the tangential stresses vanish on the fluid 
cell surface. 

The continuity and momentum equations for a variable 
viscosity system in conservative term can be expressed as 

V . V = 0 (3.16) 

pV.VV = -Vp + '7(r) (VV + VV T )) (3.17) 
where p is the non-gravitational pressure, V is the velocity 
vector; p and f) are liquid density and apparent viscosity 
respectively. The constitutive equation for the rheological 
behaviour of a power law type fluid is expressed as 
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(3.18) 


jhere 


A = V + V . 
a J a , J J , a 


and 


i 0 = r r a a = a + a* +a* + 2 a 

5 rr 66 
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(3.19) 


1 J 

For the free surface cell model, the appropriate boundary 


conditions on the particle surface (r = R^ ) in the spherical polar 
coordinates are: 
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At the cell surface, r = R 


V = 0 

r 


r<9 


* ( ~ > 


(3.20c) 


(3 . 20d) 


The radius of the fluid cell, R c , and the voidage of the 
assemblage <e are related by the expression 
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(3.21) 


The variables in the above equations are non-dimensional ised as 


v * * 

~ , v = Rp v , P 


p 0. 


c to 




V) 


U 


( IT* 


<xs ^n- 


r 


* 

A. . 
1 J 


1 J 


U 

( B 2 > 


F“ 

K P 



44 


D , 

00 2 

<• IT - ' 
Kp 


* 

T = 


u 

K ( “ ) n 

P 


(3.22) 


Uaing these non-diaens ion&l ised variables, Eqns . (3.16) to (3 21) 

are rewritten in their dimensionless for® as 
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R c (i.e R c /R p ) 
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(3.27a) 

(3.27b) 

(3.27c) 
( 3 . 27d) 


The solutions of Eqns. (3.23) and (3.24) for the boundary 
conditions stated in Eqns. (3.27) represent the detailed flow 
field (u, p , v) around a particle in the assemblage. 

The dimensionless drag coefficient defined as 


C D * V[ I p " E P U «] 



may be obtained by integrating the viscous and the pressure forces 
over the sphere surface The pressure drag coefficient is 
obtained as 
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The friction drag coefficient is obtained by integral ing the 
normal and tangential components of the viscous stress tensor as- 
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The total 

drag coefficient 

is 

simply 

the 

sum of 

pressure drag 


coefficient and the friction drag coefficient, given by the 
expression 


C = C + C 
D DP DF 


(3.30) 


3.5 SOLUTION METHODOLOGIES 

An exact solution of these equations of motion has so far 
been obtained only for creeping flow of Newtonian fluids (See 
Appendix-I). However, analytical solutions are not possible for 
non-Newtonian fluids even under creeping flow condition, owing to 
the non-linear form of the rheological equation of state (Eqns. 
2.13 or 2.14). Therefore, existing analyses of the problem 
pertain to varying degrees of approximation . Two distinct 
approaches have been used in seeking approximate solutions to this 


class of problems. 
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(1) For the steady, creeping flow of incompr essible Generalized 
Newtonian fluids, variational principles have been extensively 
used This approach yields upper and lower bounds on the 
integrated quantities such as drag coefficients and does not 
provide any insight into the physical features of the flow field 
Furthermore, in tills appt oaih , eithet the equation of continuity 
(upper bound) or the moment urn equation (lower bound) is satisfied 
It has also been shown that this technique yields true bounds only 
In the cae of Newtonian and power law fluids (Slattery, 1972) An 
important feature of such analysis is that the two bounds diverge 
mci eas ingly as the extent of non-Newtonian behaviour increases 
Excellent discussions on the variational principle approach (with 
special reference to its applications to fluid mechanical 
problems) have been given by Powlowski (1954), Bird (1960), 
Johnson (1960,61) and Leonov (1988) Similar applications of this 
technique to non-Newtonian fluid-particle systems have been 
illustrated by numet ous workers (e g see Bird, 1960, Yoshloka and 
Adachi, 1971, 1973, Slattery, 1972, Uasserman and Slattery, 1964; 
Hopke and Slattery, 1970, Mohan , 1974, Chhabra et al , 1980; Mohan 

and Raghuraman, 1976a, b, Chhabra and Raghuraman , 1984, Jarzebski 

and Malinowski, 1986; Gummal am and Chhabra, 1987; Leonov, 1988; 
Prasad et al , 1990) 

(ii) The second class of solutions rests on the 1 xne&r izt ion of 
the non-linear viscous terms in the momentum equations for 
creeping flow conditions. The essence of this approach lies in 
approximating the second i nv ar i ant of the rate of deformation 
tensor by the corresponding Newtonian value, and this together 



with none other assumptions, allows a closed form solution to be 
obtained. Perhaps Hirose and Moo-Young (1969) were the first to 
employ this approach to obtain a closed form expression for drag 
coefficient of a bubble moving in a power law fluid Since then, 
it has been extended to encompass swarms of bubbles in power law 
and Bingham plastic fluids (Bhavraju et al . , 1978), single as well 
as assemblage of solid spheres moving in power law liquids 
(Acharya et al . , 1976; Kawase and Ulbrecht , 1981a, b) . Owing to 
the nature of the simplifications made in this approach, one would 
intuitively expect the results to be applicable for weakly 
non-Newtonian conditions only. 

Beyond creeping flow conditions, no approximate solutions are 
possible even for Newtonian fluid. Numerical Methods (Leclair and 
Hamlelec, 1968), perturbation analysis (El-Kaissy and Homsy , 1973) 
and boundary layer approximations (Nishimura and Ishii, 1980a, b) 
have been employed to analyse Newtonian fluid flow through 
multi-particle systems at intermediate and high Reynolds numbers. 
For non-Newtonian liquids the application of numerical methods 
such as finite difference, finite element and boundary element 
method has been limited to the creeping flow past a single sphere. 
These methods have not been extended to multi-particle systems. 
There is a paucity of full scale numerical results for free 
surface cell model in the intermediate Reynolds number regime 
particularly for non-Newtonian flow. An approximate solution was 
obtained (Kawase and Ulbrecht, 1981a; Hua and Ishii, 1981) by 
using the Happel ’ a free surface cell model in conjunction with 
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boundary layer theory for flow of power law fluid past an 
assemblage of spheres. 

In most of the available numerical studies on fluid flow 
through particle assemblages stream function-vorticity approach of 
finite difference technique has been extensively used (Leclair and 
Hamielec, 1968 ; Hamielec et al . , 1967 ; Tiefenbruck and Leal, 1982, 
Sundarara jan and Ayyaswamy , 1984 ). Over the past two decades the 
finite element method (FEM) has emerged as a very effective method 
for solving a wide range of flow problems (Gallaghar et al , 1978, 
Baker, 1985). The creeping motion of a single sphere in Newtonian 
and non-Newtonian media has been analysed with the help of the 
finite element method (Crochet, 1982; Gu and Tanner, 1985) and the 
boundary element method (Bush and Phan- t hi en , 1984). In fact, for 
simulating low Reynolds number flows, the FEM solution procedures 
offer several attractive features over the conventional numerical 
techniques. These include, easy handling of stress boundary 
conditions, Incorporation of complex viscosity variation In the 
fluid, convenient application to complex geometry and the 
development of general codes for a variety of flow problems. 
Furthermore, the incorporation of primitive variables (u, p, v) 
approach in these studies Is definitely more convenient and 
advantageous than the stream function-vorticity formulation. It 
is well known that for the solution of the vorticity transport 
equation which constitutes an elliptic boundary value problem, a 
linear combination of vorticity and its normal derivative Is 
required to be specified on a closed boundary. In general, 
vorticity is not known cl priori along the boundry . Moreover , 
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owing to the variable viscosity of non-Newtonian fluids, the 
transformation of the equations of motion into the classical 
equation for the convective diffusion of the vorticity is vt 'y 
cumbersome. Besides, the pressure r eco\ ery process from the 
stream f unct lon-vort 1 c i ty solutions is susceptible to large 
numerical errors (Crochet, 1984) as it involves the evaluation of 
higher order derivatives of stream function and their integrtion 
along the particle surface. It is, therefore, not at all 

surprising that the primitive variable approach of the finite 
element method for solving flow problems has found widespread use. 
The penalty function approach (Bercovier and Engleman, 1979) has 
also been applied with great success to the solution of 
Navi er -Stokes equations for incompressible Newtonian fluid flow 
In this approach, the pressure variable is eliminated and only 
velocity components are solved which results in considerable 
savings m computing time However, in the case of non-Newtonian 
fluid flow, the penalty approach has not been found very useful 
(Bernstein et al . , 1981; Nakazawa et al . , 1982). In this study 
therefore, the primitive variable formulation of the flow 
equations has been adopted. In the following section, construction 
of a finite element scheme based on such a formulation is 
elaborat ed . 

3. 6 FEM FORMULATION 

The finite element method is a numerical procedure for 
solving physical problems governed by a set of differential 
equations or a variational principle. It has two characteristics 
that distinguish It from the other numerical procedures. 
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(a) It utilizes an Integeral f ormulat ion to generate a system of 

algebraic equations, and 

(b) Continuous piecewise smooth interpolation functions are used 

for approximating the unknowns . 

Amongst the various mathematical theories, the variational 
principle has been (and continues to be) the most popular basis of 
forming a finite element approximation. However, for many 
problems in fluid mechanics, variational principles do not exist. 
This la so because in the Euler ian frame of reference, the 
expression for conservation of momentum is explicitly non-linear 
In such cases, the weighted residual approach offers a more 
powerful tool for constructing the integral formulation Here, 

the residue of each governing conservation equation is itself not 
set equal to zero every where but its integeral with respect to a 
selected weighting function is required to vanish (Finlayson, 
1972). Several approximation techniques based on this approach 
are available such as the collocation method, the Galerkin method, 
the least square method etc. The Galerkin's weighted reeidual 
technique is by far the most versatile and popular method (Baker, 
1973; Popinski and Baker, 1976; Barrett and Elliott, 1982; 
Engleman and Sanl , 1986) and is also more suitable for elliptic 
operator problems. 

In present work, the Galerkin’s weighted residual approach 
has been used to construct a finite element scheme to solve the 
governing equations. At the outset. It is worthwhile to recognize 
that owing to the axisymraetry in 4 > direction, it Is much more 
convenient to numerically solve the equations of motion along with 
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Fig 3 2 Typical spatial discretization of flow domain (« = 0 999) 
utilising quadrilateral elements with exponential step 
size in radial direction 


(U5.P5.V5) 


3 



o X 


Parabolic eight noded element showing typical node 
numbering and local coordinate system. 
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the associated boundary conditions (which have been presented in 
Section 3 3 in spherical poloar coordinates) in a cylindrical 
coordinate frame (y , , x) shown in Fig 3 1 Indeed , the FEW 

formulation is considerably simplified due to the aeixymmetry in ^ 
direction and has a form which is only slightly different from a 
2-D formulation. 

The spatial discretization of the flow domain (Fig. 3.2) has 

been done using 8-noded quadrilateral elements (Fig 3.3) of the 

serendipity type containing one node at each corner and one at the 

midpoint of each side. The variation of V , V and p over each 

x y 

element is defined as 
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(3.31) 


where V , , V . and p. are the nodal values of the variables. Also 

and are the int erpolat ion functions (shape functions) for 

velocity and pressure respectively. Since V and V have a 

x y 

similar range of variation and satisfy similar equations, same 

shape functions have been used for both V and V . One of the 

x y 

problems usually encountered in solving the Navi er-Stokes 
equations by the finite element method has been the occurrence of 
spurious pressure fluctuations for certain choices of elements 
(Taylor and Hood, 1974; Olson and Tuann , 1976; Jackson and Cliffe, 
1980; Sani , 1980). The "mixed interpolation” method obviates this 
difficulty, though one can not say that the problem has been fully 
understood (Rae, 1982). It is now a standard practice (Taylor and 
Hughes, 1984) to approximate the variation of pressure by shape 


functions of one order lower than those used for velocity. 


The 
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pressure In the present scheme is interpolated by means of a 
bilinear polynomial, whereas the velocity components are 
approximated by serendipity shape functions giving an 

approximation of P^-C° type These features are implemented by 
using the same element geometry, where all the eight nodes are 
associated with velocities but only corner nodes Involve pressure. 
The various options available for the choice of the elemental 
shapes and the shape functions have been adequately discussed by 
Crochet et al. (1984). 

Dropping the asterisk for convenience from the 

non-dimensional eqns . (3.23) to (3.27), the weighted residual 

statement for the governing equations of continuity and momentum 
in X and Y directions , may be written as 

III H. (7.V) y dy dx d^, i = 1, ..., m (3.32) 

a 

JJJ N i e x- [ v - W + Vp - ^1— T) (VV + VV T ) }] y dy dx 

ft r 

= [0] , i = 1, ... 1 (3.33) 

JJJ “i V ( 

ft 

= (03, 1 = 1, ... 1 (3.34) 

where ft refers to the solution domain, m is the total number of 

pressure unknowns and 1 is the total number of velocity unknowns. 

It is interesting to note that even though the pressure gradient 
forms a part of the momentum equations, yet there is no obvious 
equation for obtaining the pressure variables. Therefore, 

pressure is evaluated only indirectly through the continuity 
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equation in the sense that when the correct pressure field is 
substituted into the momentum equations, the resulting velocity 
field satisfies the continuity equation One approach to avoid 
spurious oscillations due to the non-satisfaction of continuity 
equation is to snake the pressure field explicitly compatible with 
the continuity equation 

The choice of the shape functions and as weighting 

functions for minimising the residue of the continuity and 
momentum equations respectively, can be shown to yield a 
diagonally dominant system of matrix equations in 
viscosity-dominated situations as encountered in the applications 
involving particle assemblages. However, for convection-dominated 
situations, upwind weighting functions (in Petrov-Gal erkin method) 
which give more weightage for residue in the downstream direction 
as compared to that in the upstream direction have also been 
proposed. In the present work, since Reynolds number is not very 
high (Re < 100) the above choice of weighting functions is 

considered to be adequate and justified. 

Invoking the Green’s theorem, the second order terms in the 
momentum equations can be reduced leading to what Is commonly 
called as the weak form of weighted residual statement, as given 
below: 

m Nj (V. <7V + Vp) + V N i - gj- (W + VV T )> y dy dx dj> 

a 1 * r 

- # n . N. (VV + VV T )> y dl &4> = [0] (3.35) 

c r 

where n is the outward normal at a point on the boundary surface c 

Substituting the approximation for V , V and p 

x y 


of the domain ft 
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from Eqn. (3.31) in Eqns . (3.32) to (3.35), the final form of 
weighted residual statements for continuity equation and momentum 
equations turn out to be 
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In above equations, the terms F , and F . on R.H.S. are the 

x , i. y , i 

force contributions from the deviatoric stress on the cell 


boundary . However, in the present work, there is no need to 
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evaluate these forces since, a simpler method has been adopted 
(which is described subsequently) for incorporating the cell 
boundary conditions of the free surface cell model 

Equations (3.37) and (3.38) exhibit strong non-linear 
behaviour both due to the presence of the inertial ter#a and the 
shear dependent viscosity . To facilitate the ease of obtaining a 
solution , these non-linear terms have been quasi 1 in ear ised , by 
evaluating them on the basis of the solution corresponding to the 
previous iteration. For the sake of clarity , the use of the 
previous iteration solution is indi cat ed with an overbar (e.g. V . 
and V wherever it occurs explicitly The elemental integrals 
in the Equations (3.36) to (3.38) have been numerically evaluated 
using Gauss Legendre quadrature with 3x3 sample points Thus, 
the system of the non-linear int egr al equations is reduced to a 
system of non-linear algebraic equations . In matrix form, for an 
element , these equations can be expressed as 

[K] { 4 > x ) = (F> (3.39) 
where [K] is a 20 x 20 coefficient matrix, } is the vector 
containing al 1 the nodal unknowns e.g u^ , p^ and v^ and { F ) is 
the right hand side vector also known as flux vector . Non-zero 
values of the components of flux vector arise only f or elements 
lying on the boundary as the contributions for the interior 
elements mutually cancel each other. 

Before discussing the assemblage of these elemental 
coefficient matrices into a Global Coefficient matrix, the 
evaluation of the non-Newtonian viscosity and the method of 
impl ement ing the free surface cell boundary conditions is 
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presented In the following section. 

3. 7 Evaluation of Non-Newtonian Viscosity 

The non-dimensional viscosity for a power law fluid is given 
by 

n-1 



where I^, the second invariant of the rate of deformation tensor 
in cylindrical coordinate system (y, 4 > , x) with axisymmetry 
prevailing m 4 > direction, can be expressed as: 
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(3.40) 


The viscosity is determined at the Gauss points for evaluating the 

integrals in equations (3.36) and (3.37) for any interior element 

The velocities V and V as well as the global coordinates x and y 

x y 

at the Gauss points are determined using the interpolation 
formulae 
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Here V . and V . represent the velocity components at the nodal 
xk yk 

points obtained from the previous iteration. While determining 
the shear stress and the viscous drag on the particle surface, the 
viscosity is evaluated at three gauss points located on the side 
of element which lies on the particle surface. 

For the determination of the second invariant (I 2 ) and also 
other terms contained in equations (3.36) to (3.38), the first 
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or aer der i vat i ves of shape functions with respect to the global 
coordinates (x , y ) are needed . These are obtained as follows : 
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(3 43) 


where J is referred to as the Jacobian of transformation and £ ^ , 
are the local curvilinear coordinates of an element. The 
Jacobian matrix can be evaluated explicitly by substituting the 
variations of x and y from equation (3.41b) into the definition of 
[ J ] . Thus , 
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(3.44) 



60 


3. S Impl*«»nlation of Boundary Conditional 

The non-dimensional boundary conditions for the free surface 
cell model in cylindrical coordinat e system are written (without 
asterisk) as 
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where V is the normal velocity and r . is the shear stress on the 

cell surface It is appropriate to add here that while applying 

the stress boundary condition at the cell surface, it is 

convenient to apply it in polar spherical coordinate system. Thus 

only and V & have been transformed into cylindrical velocities, 

V and V while applying the cell boundary conditions . Since one 
x y 

of the local coordinates (£ ^ or £ 2 depending on respective side of 

the element) coincides with the radial direction, the first order 

derivative in this direction can easily be converted into first 

order derivative with respect to the local coordinate system 

Substituting the elementwise interpolation of V and V for 

x y 

each element lying adjacent to the cell boundary, the zero shear 
stress boundary condition i.e. Eqn (3 45c) becomes 



(3.46a) 
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where { ^ has been taken to lie in the radial direction (r ) , with 
no loss of c ener&l ity The zero normal velocity at the cell 

boundary can be written in the form 

V cos & + V sm 0 = 0 (3 46b) 

xj yj 

where j is the number of each node which lies on the cell 

boundary. For nodes lying on the cell boundary. Equations (3.46a) 

and ( 3 46b) replace the respective momentum equations at these 
nodes For the nodes on the particle surface the essential 

boundary conditions of the form 

V . = 1.0, V , = 0.0 (3.47) 

xj yj 

are applied. 

3.0 GLOBAL MATRIX AMD SOLUTI OM TECHNIQUE 

After the formation of an elemental coefficient matrix [K] 
and the associated right hand side vector { F } in Eqn . (3.39), 

further logical steps in arriving at the final solution are: 

(a) Assemblage of el emental matrix into a Global matrix . 

(b) Suitable matrix operations to obtain the detailed flow field 
solution. 

The method adopted for solving the assembled matrix equations 
has a significant bearing on the computer storage requirement and 
execution time. Uhen a large set of equations are to be solved, a 
simple matrix inversion technique will be very cumbersome and 
expensive due to large computer core memory requirement. In such 
conditions the direct elimination FRONTAL SOLUTION METHOD has 
proved to be very ef f ect ive for solving positive definite 
matrices. The size of the G 1 obal mat r ix in present problem is of 
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the order of 1000 x 1000 and therefore, the choice of FRONTAL 
method as a matrix solver becomes inevitable. 

The FRONTAL method, originally developed by Irons (1970) has 
been based on direct Gaussian elimination for solving symmetric 
matrices where the leading diagonal is always used as the pivot. 
For unsymmetric matrices which are encountered in a wide range of 
problems in engineering fluid mechanics, the most suitable pivot 
is not necessarily on the leading diagonal and FRONTAL solution 
routines exist for off diagonal pivoting (Hood, 1976). These, 
however, tend to be more time consuming and therefore, for the 
solution of non- symmetric matrices, in the present work only 
diagonal pivoting has been employed. The philosophy and operation 
of a FRONTAL method has been covered at length by Taylor and 
Hughes (1981). The method excellently manages the core memory 
requirement by writing the fully assembled nodal equations on to 
the disc memory as soon as contributions from all the elements 
which surround that node are added. The assembled equations are 
reduced by Gaussian elimination before being written on tape. 
After all equations have been completed, the back substitution 
process of unknowns is performed by recalling the equations in 
reverse order from tape. Thus, since assembled equations are 
never stored in core memory, the FRONTAL method drastically 
reduces the size of memory to perform all the required matrix 
operations . 

The non-linearities of the Eqns . (3 37) and (3.38) earlier 
mentioned necessitate an iterative type of solution procedure. It 
was found that under some circumstances, the solutions have to be 
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under-relaxed between two successive iterations to increase the 
stability of the numerical results . The computational work 
reported herein was performed on DEC 1090 and ND 560 systems The 
CPU time taken for a typical iterat ior using 18 x 8 Mesh varied 
from 3 to 4 minutes 

3.10 POST PROCESSING OF RESULTS 


After obtaining converged flow field solutions, the results 
are processed further to evaluate various quantities of interest 
such as the dimensionless stresses at the particle surface, the 
drag coefficients (C Dp , C Q p and C D ) and the dimensionless stream 
function value at each node. For the sake of convenience in 
presentation, two sets of dimensionless stresses are defined using 
the expressions 


T 


r Jp\5 2 

oo 


P 


P/PU 2 

0D 


(3 48a) 


and , t ' = r / ( 


i?U 


m 


i?U 


); 


p/C 


00 


(3.48b) 


P P 

The use of t and p is advantageous for power law fluids, 
since the scaling quantity used in their def inition Is independent 
of the power law index n. On the other hand, t ’ and p ’ are useful 
for Newtonian fluids with constant viscosity, since both t' and p’ 
will be of the order of 1 for the range of Reynolds numbers 
considered in the present study. In terms of the numerically 


* * * 


predicted flow field quantities V , V and p , the dimensionless 

y 

stresses are given as - 



P 


(3.49a) 



64 


r 

Ie“ 


if" { * 

f' ^ 


. ,v s 


ar 


} 


r = l 


(3.49b) 


«?TV 


The velocity gradient 


represented as 


3 


in eqn . (3.49b) can , in turn, be 


ar 


av 


Sr 


Ar 

1 


1_ S_ 


* * 

( - V sxn & + V cos & ) 
x y 


1 

8 <9N 


* r Tsv-i ( - V . Sin © + V . cos ©) 
Ar j = i *■> yj 


(3.50) 


where, Ar is the dimensionless radial step and £ ^ is the local 

coordinate in the radial direction. The summation in Eqn. (3.50) 

is carried over all the nodes of the element of which the surface 

node under consideration forme a part. 

The drag coefficients and can be evaluated by 

DP DF 

integrating equations (3.28) and (3.29), via 3-point 
Gauss -Legendre quadrature method along the boundary curves of 
elements which lie on the particle surface. Thus, eqns. (3.28 and 
3.29) can be rewritten as: 

e=E 


'DP 


E J E «i Pj. sin 


6 cos 3 6.3 


(3.51a) 


e=l D e i-1 


e=E 


'DF 


e “ c s 8 SN * 

K E J E * { 2 -4 V 

Ke e=l ft e i = l 1 ax 


S N 


cos 3 + 


S Y 


i * 

=- V . sm 3 
* xi 


SU 


ax 


t V . sin 3 
* yi 


} ein 


3 63 


(3.51b) 


where the subscript ft signifies integration over the boundary of 
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a typical surface element O e and the summation is carried out for 

all the elements lying adjacent to the particle , with E denoting 

8 

the total number of such elements. 

For the sake of plotting the streamlines, dimensionless 
Stream function (y ) values have been calculated at each node by 
integrating the expression 


V 

r 



W 


(3 52 ) 


with respect to angle & at a constant r The stream-function on 

the axis (& = 0) is assigned the value zero and marching is done 

in the & direction from B = 0 to & = n m steps of LB at a 

constant r, using the FEM grid for the problem. This procedure is 

* 

repeated for each value of r ranging from r = 1 to r = R in steps 
of Ar , to obtain the stream function values for the whole domain. 



CHAPTER IV 


VALIDATION OF NUMERICAL SCHEME 

In order to establish the correctness of the numerical scheme 
employed in the present work, the drag coefficients and the 
detailed flow field solutions have been compared with existing 
analytical and numerical results in literature. In addition to 
these comparative studies, the effects of various numerical 
parameters, namely, the radial step size, the grid pattern and the 
initial guess profiles, etc., on the accuracy and rate of 
convergence of the predicted results have also been investigated. 

4.1 NEWTONIAN FLUID FLOW PAST MULTI -PARTICLE ASSEMBLAGES 

Numerical results using FEM have been obtained for Newtonian 
as well as non-Newtonian fluid flows relative to multi-particle 
assemblages at low and intermediate Reynolds number regimes. 
However, a detailed bench mark comparison is possible only for the 
creeping flow of a Newtonian fluid, for which an analytical 
solution using Stoke’s stream function was earlier provided by 
Happel (1958). The drag coefficients evaluated for any other flow 
situation using Happel cell model is based on certain 
approximations and linearization techniques. In Table 4.1, a 
comparison between the analytical results reported by Happel 
(1958) and the present numerical results for the creeping flow of 
Newtonian fluids past an assemblage of particles is presented. 
Though the analytical expression for the total drag coefficient 
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hass been explicitly provided by Happel , the correspc ding 
expressions for friction and pressure drag coefficients are, 
however, not available in open 1 iterature and hence are Included 


in Appendix 

I An examination 

of Table 

4.1 

shows that 

the 

individual 

contributions (i e. 

friction 

and 

pressure 

drag 


coefficients) as well as the values of total drag coefficients 
obtained herein are in complete agreement with the analytical 
results over the entire range of values of voidage the maximum 
deviation observed between the two values is of the order of 1.3% 
for the highest value of voidage (« = 0.9999) used here In view 
of this, it is therefore safe to conclude that the correspondence 
between the analytical results and present numerical results is 
excellent . 

In Tables 4 2 and 4.3, a comparison between the analytically 
and numerically obtained surface pressure and shear stress values 
for creeping flow are presented at different angular positions for 
two voidages (e = 0.4 and 0.9). These two voidages typically 
represent concentrated (« = 0.4) and dilute (e = 0.9) assemblage 
of particles. It is evident that there is an excellent agreement 
between the analytical and numerical predictions of surface shear 
stress at all angular positions. For surface pressure, it is 
observed that although there is almost total agreement between the 
numerical and analytical results (Maximum deviation < 1%) over a 
wide range of angular positions (20° < 9 < 180°), there exists a 
narrow zone in the'front region (0° < 9 < 20°) over which the two 
results show a higher deviation (maximum about 10%). The reason 
for this could be that the pressure in this narrow zone is much 
smaller than the highest pressure occurring in the computational 
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Table 4.2 Compar Ison o £ Analytical and Numerical Results 

- Non dimensional surface pressure. P / (yjU /R ) and 

& CD p 

shear stress t /(r;U /R ). 

S O) p 

« = 0 4 , Re = 0 001 n = 1 


Angular 

posl t ion 


P 

s 

T 

S 


Analytical 

Numerical 

Analytical 

Numerical 

0 

0 . 0 

0 0 

0 0 

0 0 

10 

-0 8199X10 4 

-0 888x1 0 4 

-0 198x10 s 

0.200xl0 S 

20 

-0 3255xl0 S 

-0 . 331xlO S 

0.390xl0 5 

0.391xl0 S 

30 

-0.723xl0 5 

-0.728xl0 5 

0.570xl0 5 

0.568xl0 S 

40 

-0 1262xlQ 6 

-0. 127xl0 6 

0 733xl0 5 

0.730xl0 5 

50 

-0 1928xl0 6 

-0.1 93xl0 6 

0.874xl0 S 

0.870xl0 S 

60 

-0 2699xl0 6 

-0.27 0x10^ 

0.988xl0 S 

0 984xl0 5 

70 

-0 3551xl0 6 

-0.355xl0 6 

0.107xl0 6 

0. 106xl0 6 

80 

-0.446X10 6 

-0.446X10 6 

0.1124xl0 6 

O.lllxlO 6 

90 

-0.5398xl0 6 

-0.539xl0 6 

0. 1141xl0 6 

0.113xl0 6 

100 

-0.6335xl0 6 

-0. 633xl0 6 

0 .1124xl0 6 

0 UlSxlO 6 

110 

-0. 7244xl0 6 

-0.724xl0 6 

0.1072xl0 6 

0. 1068xl0 6 

120 

-0.809 7x10^ 

-0.809xl0 6 

0.9885xl0 S 

0.984xl0 S 

130 

-0.8868xl0 6 

-0.886xl0 6 

0 8743xlO S 

0.870xlO S 

140 

-0.9533xl0 6 

-0. 952xl0 6 

0 7337xl0 5 

0.730xl0 S 

150 

-0. 1007xl0 7 

-0. lOixxlO 7 

O.57u7xl0 S 

0.568xl0 S 

160 

-0. 1047xl0 7 

-0.105xl0 7 

0.3903xl0 S 

0.390xl0 S 

170 

-0.107 lxl 0 7 

-1 107xl0 7 

0.1982xl0 5 

0.200xl0 S 

180 

-0.1079xlU 7 

-0 . 108x1 0 7 

0.0 

0.0 
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Table 4 3: 

Cofnpar Ison 

of Analytical 

and Numerical 

Results 


- Non dimensional surface 

pressure, P 

s 

and 


shear f tress r / (y> U /R ) 

S top 

& = 0.9; Re = 0.001 si = 

i 


Angi lar 

position 

9° 


p 

s 

T 

& 


Analytical 

Numerical 

Anr lytical 

Numerical 

0 

0.0 

0 0 

0 0 

0 0 

10 

-0 1520x10* 

-0 .139x10* 

0.1565x10* 

0.1564x10* 

HO 

-0 6035xlQ 4 

-0.588x10* 

0.3082x10* 

0.3080x10* 

30 

-0 1340x10* 

-0.132x10* 

0.4506x10* 

0.4500x10* 

40 

-0.2341x10* 

-0 232x10* 

0.5793x10* 

0.5780x10* 

50 

-0.3575x10* 

-0 356x10* 

0.6904x10* 

0 .6880x10* 

60 

-0.5004x10* 

-0 499x10* 

0.7805x10* 

0 .7780x10* 

70 

-0.6585x10* 

-0 657x10* 

0.8469x10* 

0.8440x10* 

SO 

-0.8270x10* 

-0.825x10* 

0.8876x10* 

0.8860x10* 

90 

-O.lOOOxlO 5 

-0 999x10* 

0.9013x10* 

0 8980x10* 

100 

-0.117 4x1 0 S 

-0.117xl0 S 

0 8876x10* 

0.8860x10* 

110 

-0.1342xl0 5 

-0.134xl0 S 

0.8469x10* 

0.8440x10* 

120 

-0.1501X10 5 

-O.lSOxlO 5 

0 7805x10* 

0 7780x10* 

130 

-0.1643xl0 S 

-0.164xl0 S 

0.6904x10* 

0 6880x10* 

140 

-0.1767xl0 S 

-0.177xl0 S 

0.5793x10* 

0.5780x10* 

150 

-0.1867xl0 S 

-0.187xl0 S 

0.4506x10* 

0.4500x10* 

160 

-0.1941xlO S 

0 194xl0 S 

0.3082x10* 

0.3080x10* 

170 

-0.1986xl0 S 

-0 198xlO S 

0.1565x10* 

0.1564x10* 

ISO 

-0.2001xl0 S 

-0 2000xl0 S 

o 

o 

o 

o 
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domain. It is a well known prooleo to the numerical analyst that 
during the simultaneous solution of the (u, p, v) variables, the 
smallest values of pressure undergo wild fluctuations. In order 
to obviat e this difficulty of spurious pressure oscillations , one 
of the commonly adopted procedures is to exempt the lower values 
of pressure from the convergence test. In the present work also, 
a cut-off value for pressure was prescribed below which the 
convergence test was not applied It is believed that the higher 
deviation between the analytical and numerical results for 0 i & 
£ 20°, arises from the occurrence of pressure oscillations m that 
region. However, the pressure drag coefficient remains uneffected 
by these oscillations because it is an integrated quantity on 
which the contribution due to very low values of pressure is 
insignificantly small. Notwithstanding these aspects, it is 
readily recognized that the analytical and numerical values of 
surface stresses are in excellent overall agreement 

The validity of the numerical results has been finally 
ascertained by comparing the numerically computed radial velocity 
and angular velocity values against the corresponding analytical 
values in Tables 4.4 and 4.5 Table 4 4 corresponds to a dense 
assemblage of particles (€ - 0.4), whereas Table 4.5 contains the 
results for a dilute assemblage (e = 0.9). The velocity values 
have been compared at several radial positions Excellent 
agreement between the analytical and numerical results is evident 
from the Tables. 

In view of the above discussed comparisons between the 
present results and the analytical results of Happel (1958), the 
credibility of the numerical scheme adopted and the validity of 



Comparison of analytical and numerical results - non-dimensional 
radial velocity and angular velocity. 
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Table 4 5. Comparison of Analytical and Numerical Results - Non-dimensional 
radial velocity and angular velocity. 

« = 0.9, Re = 0.001, n = 1.0 


73 






lO 

co 

oo 

UO 


rH 

03 

m 

HD 

* 




off 

CO 

O' 

Os 

Os 

03 

CO 

Os 

o 

u 

1 


o 

Os 

o 

CO 

«cf 

CD 

Os 

HD 

o 

tH 


1 

<3> 

o 

sh 

vO 

03 

rH 

O 

rH 

CO 

m 

so 

6 

1 

> 

• 

. 

. 

. 

• 


. 

. 

» 

• 


! 

j 


rH 

o 

o 

o 

0 

1 

o 

o 

o 

o 

a 

mm 

o 

1 


1 

1 

i 

s 

1 






o 

1 












Os 

1 



vO 

03 

03 

O' 

03 

o 

VO 


tn 


1 



*S3> 

O' 

co 

CO 

O 

03 

CO 

CO 

Os 

• I! 

1 


o 

Os 

O 

CO 

K? 

Os 

Os 

HD 

rH 

tH 

pH 

1 

<b 

o 


vO 

03 

rH 

o 

rH 

CO 

tn 

V© 


1 

> 

• 

. 

♦ 

. 


. 


• 


• 

C 

1 


tH 

o 

o 

o 

O 

o 

c 

o 

o 

o 

< 

! 

1 


1 

1 

1 

1 

1 









o 

■H 

03 

CO 

*3 

to? 

"C? 

Os 




1 


vO 

CD 

O' 

Os 

Os 

O' 

vO 

oo 

03 

V© 


1 


VO 

CO 

03 

tH 

03 

O' 

so 

rH 

tH 

V© 


1 

■ft 

CO 

vO 

in 

03 

•pH 

o 

tH 

CO 

*3* 

CO 

pH 

1 

> 

• 

. 



* 




«o* 

un 

<d 

1 


O 

o 

o 

o 

Q 

o 

o 

O 

• 


u 

1 


! 

! 

i 

1 

1 




O 

o 

-H 

I 












u 

1 












4) 

1 



03 

CD 

o 

O 

rH 

CD 

VO 

o 

tH 

6 

1 


O 

if) 

tH 

VO 

so 

Ch 

Os 

03 

tn 

o 


1 


O 

CS 

co 

CO 

O 

CO 

Os 

03 

*3* 

H 

52 

1 

U 

if 

<3 

■SO® 

-5? 

«5f 

co 

03 

03 

tH 

O 


1 

> 

• 


• 


* 


» 

• 

• 


O 

1 


e 

O 

O 

o 

o 

O 

o 

o 

o 

o 

o 

1 












o 

1 







CO 






1 


o 

rH 

CO 

Lf! 

CD 

tH 

CO 

03 

HD 

in 

n 

1 


vO 

CD 

lD 

Os 

CO 

00 

vO 

Os 

«a* 

HD 


1 


SO 

CD 

03 


03 

H 

HD 

tH 

•a* 

CO 


1 


CD 

HD 

tn 

03 

tH 

o 

rH 

CO 


in 

rH 

1 

> 

• 

. 


. 

• 

• 

• 



• 


1 


O 

O 

o 

O 

Q 

o 

o 

O 

o 

o 


1 


1 

1 

1 

1 

1 






*H 

! 












-H 

1 




vO 

Os 

CO 

CO 

tH- 

03 

o 



l 


o 

«Cf 

▼H 

If 

lO 

SO 

Os 

03 

tn 

o 

pH 

1 


o 

Os 

00 

CO 

o 

CO 

03 

03 


o. 

<d 

1 

u 

lO 


"53* 

<5? 


co 

03 

03 

tH 

o 

c 

1 

> 

» 

. 



• 

• 

* 

» 

* 

• 

< 

1 

1 


o 

o 

O 

o 

o 

o 

O 

O 

O 

o 





03 

«ctf 

Os 

Os 

Os 

rH 

rH 

CO 

00 


1 


o 

O' 


CO 

<3* 

*3“ 

vO 

<«3> 

<=? 

Os 


l 


o 

o. 

O 


O'. 

«3* 

Os 

CO 

if 

o 


1 

<h 

lO 

CO 

CO 

tH 

o 

o 

o 

tH 

03 

co 

pH 

1 

> 

» 

t 

• 


■ 

« 

• 

• 

* 

* 


i 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 


! 


1 

1 

l 

1 

1 






*H 

\ 













1 













1 


o 

o* 

if] 

tH 

03 

Os 

CO 


03 

tn 

e 

I 


vO 

o. 


if 

00 

CO 

Os 

ud 

rH 

•H 

3 

\ 

u 

SO 

m 

CO 

if 

o 

CO 

rH 

co 

tn 

03 

52 

S 

> 

CD 

CO 

CO 

O' 

tH- 

if 

if 

CO 

03 

rH 

0 

1 


• 

• 

• 

• 

• 


• 

• 


* 

o 

1 


O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

00 

l 












II 

! 

1 



on 

VO 

■H 

CO 

rH 

o 

co 

tH 

IH 

pH 

! 


o 

O' 

CO 


■O' 

if 

vO 

*5? 

HD 

Os 


i 


o 

Os 

o 


Os. 

< 53 * 

Os 

CD 

in 

o 

u 

1 

> 

If] 

C0 

CO 

tH 

o 

o 

o 

rH 

03 

CO 

►H 

1 


* 

. 

* 

• 

* 

* 

* 

* 

• 

* 


1 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

>> 

\ 


1 

{ 

1 

1 

1 






pH 

1 













1 












CJ 

1 


o 

03 

03 

tH 

Os 

«ar 

03 

CO 

tH 

03 

< 

1 

u 

vO 

O' 

*cf 

m 

03 

co 

aD 

«a* 

rH 

rH 


1 

> 

sO 

if) 

CO 

f 

o 

CO 

rH 

oo 

if 

03 


1 


CO 

CD 

co 

r- 

Os 

tn 

tn 

co 

03 

*rH 


1 


o 

O 

o 

o 

o 

o 

o 

o 

O 

o 





Os 

tH 

rH 

tH 

Os 

O' 

m 

O' 

vO 




o 


O 

tH 

£H 

Os 

vO 

rH 

O' 

tn 


U 


o 

o 

tH 

03 

03 

CO 

**& 

HD 

O'. 

Os 




tH 

rH 

tH 

tH 

tH 

tH 

rH 

tH 

tH 

rH 


03 


153 0 00 0.3504 0.00 0.3504 0.00 0 6069 0 000 0 6070 0.7008 0 7008 



74 


results obtained from it, are firmly established atleast for 
creeping Newtonian fluid flow. In contrast to this, the validity 
of the numerical results for the intermediate Reynolds number 
regime cannot be made in an equally firm manner . The reason for 
this is that an exact solution f or flow past a multi-particle 
assemblage has not been obtained till date due to the 
non-linearity associated with the inertial terms . Furthermore , 
very little theoretical work has been reported on the modelling of 
fluid flow and heat and mass transport processes occurring in 
mult i -particle assemblages outside the creeping flow conditions 
using the Happel’s free surface cell model. El-Kaissy and Homsy 
( 1973 ) extended Happel ’ s free surface cell model beyond creeping 
flow regime by carrying out a regular perturbation analysis . They 
provided analytical expressions for the friction factor as a 
function of the bed voidage (0.3 < « < 0.6) and Reynolds number . 
The comparison of present numerical results with those of 
El-Kaissy and Homsy is shown in Table 4.6. The relation between 
the total drag coef f ici ent and the friction factor f or flow 
through a fixed bed may be establ ished in following way. The 
Blake-Kozeny friction factor is expressed as 


* 3 d 
r _ Ap « p 

bk " T (l-«) ..2 

<sb 

Furthermore , by a force balance on the cell, 

Ap „ F P 

L 4/3 n R 3 

c 


(4.1) 

it can be shown that 

(4.2) 


Therefore, 


f 


bk 


4/3 n R 


3 


€ 

(l-«) 



(4.3) 



Table 4.6 Comparison of present numerical results with 
the analysis c £ El-Kaissy and Homsy, and the 
predictions o * Ergun f s equation. 


4i 

Re 


f bk 

(present ) 

f bk 

(El-Kaissy & Homsy) 

f bk 

(Ergun) 

0 . 3 

1.0 


112 289 

111 336 

106 . 755 

0.4 

1 . 0 


98 30 

98 0507 

91.75 

0 . 5 

1 0 


85 47 

85.3051 

76 50 

0 . 3 

5.0 


22 4295 

22 2676 

22.75 

0 . 4 

5 0 


19.66 

19 . 623 

19.75 

0 . 5 

5 0 


17 .10 

17.08 

16 75 

0 3 

10 0 


11 . 2350 

11 . 134 

12 25 

0.4 

10.0 


9 83 

9.832 

10 75 

0 . 5 

10 0 


8.558 

8.581 

9 .25 

0.3 

20.0 


5 . 60 

5.568 

7.000 

0 4 

20.0 


4 927 

4.9567 

6.250 

0 . 5 

20.0 


4 2944 

4 .3681 

5 500 

0.3 

50.0 


2 . 2567 

2 2307 

3.85 

0 . 4 

50.0 


1 99 

2.0968 

3.55 

0.5 

50.0 


1.7563 

1 9638 

3.25 

0 . 3 

100.0 


1.1487 

1 .1213 

2.80 

0.4 

100.0 


1 . 0363 

1 2521 

2.67 

0.5 

100.0 


0.9260 

1 .3681 

2.50 

f bk 

= 

3 

4 

_ 3 

C d * 




(present ) 
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Expressirg Drag force in terms of drag coefficient and on 
further simplification, equation (4 3) reduces to 

£ bk 1 1 c d « 3 < 4 - 4 > 

An examination of Table 4.6 cl early shows that good 

correspondence exists between the two sets of results upto Re = 

50; the maximum divergence between the two sets of results is of 

the order of 11% for the case of « = 0.5 and Re = 50, whereas the 

two values differ only by 2% or so for the case of <s = 0.3 and Re 

= 100. This trend is consistent with the findings of El-Ka say 

and Homsy (1973) who asserted that the upper limit of Re upto 

which their analysis Is applicable decreases with the Increasing 

value of voidage. 

Other theoretical results on the flow of Newtonian fluids 

through packed beds at intermediate Reynolds number which are 

based on the Happel's free surface cell model, have been due to 

Nishimura and Ishii (1980a, b). Using boundary layer assumptions, 

they reported extensive results on friction factor and Sherwood 

number as functions of bed voidage and Reynolds number and claimed 

5 

their results to be applicable in the range, 10 ^ Re < 10 . 

Intuitively, it is highly unlikely that the boundary layer flow 

approximation will be valid at values of Reynolds number as low as 
10. Several unusual features of their results have already been 
discussed In Chapter II. It is for these reasons that a 
comparison of present results with those of Nishimura and Ishii 

(1980a, b) has not been made. 



Table 4.7 Comparison of variational and present numerical results on drag 
coefficients for power law fluids (Re = 0.001). 
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4.2 NON— NEWTQMI AN FLUID FLOW PAST MULTI -PARTICLE ASSEMBLAGES 

The utility of the present numerical scheme in dealing with 
the creeping flow of power law fluids past multi-particle 
assemblage has also been examined. In Table 4.7, the present 
results are compared with the approximate upper and lower bounds 
on drag coefficient reported by Mohan and Raghuraraan (1976a). It 
is seen that a reasonable correspondence exists between the two 
sets of results obtained using entirely different techniques It 
is generally recognized that the divergence between the upper and 
lower bounds increases as the degree of non-Newtonian behaviour 
increases, i . e . , as the value of flow behaviour Index (n) 
decreases. Besides, the variational approach is tantamount to 
either satisfying the equation of continuity (upper bound) or the 
momentum equation (lower bound) only. Notwithstanding these 
aspects of the variational analysis, the agreement between the two 
sets of results Is regarded to be satisfactory and acceptable, 
albeit the present results are believed to be more accurate 
especially In the lower range of n. 

To the best of our knowledge, there have been only two 
theoretical studies relating to the flow of power law liquids 
through particle assemblages outside the creeping flow regime 
(Reynolds number greater than unity) where, the free-surface cell 
model has been employed to represent the inter-particle 
interactions. The boundary layer approximation has been used by 
Hua and Ishii (1981) and Kawase and Ulbrecht (1981b) to solve the 
governing equations for the high Reynolds number flow of power law 
liquids. The values of drag coefficients have been obtained in 
the Reynolds number range of 10 to 1000. The validity of the 
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boundary layer approximation appears to be questionable in this 
flow configuration especially in dense assemblages Moreover, the 
theoretical r es _lts reported by Hua and Ishii (1981) do not reduce 
to the expected behaviour for a single sphere as the bed voidage 
approaches unity It is to be mentioned at this stage that in the 
present study the upper reynolds number has been confined to the 
value of 20. Therefore, the numerical results were not compared 
with the approximate results of Hua and Ishii (1981) or Kawase and 
Ulbrecht (1981b). 

4.. 3 EFFECT OF NUMERICAL PARAMETERS ON COMPUTED RESULTS 

The effects of various numerical parameters, namely, the 
radial step size, grid pattern and initial guess values on the 
computed values of drag coefficient and the rate of convergence of 
the solution have been investigated. These are described in the 
ensuing sections . 

(a) Radial Step Size: 

In Table 4.8, drag results are presented for two different 
meshes which consist of 18 elements m the angular direction and 8 
or 4 elements in the radial direction. The grid spacing is 
uniform in both radial and angular directions for these meshes. 
The values of drag coefficient obtained using these grids are 
compared with the analytical results under creeping flow 
conditions, for a range of values of bed voidage. It Is readily 
seen from the results in Table 4.8 that upto about <e < 0.5, the 
two grids give virtually identical results. Only for the values 
of voidage greater than 0.5, the use of a finer grid such as 18 x 
8 mesh appears necessary. Further refinement of the mesh using 18 
x 16 elements resulted in virtually no change in the values of 
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Table 4 9 


Performance of uniform grid vs exponential grid at 
different voidages (Re = 0.00). 



n 

Cp (available 
in literature) 

Uni form gr id 
(18x8) 

Exponential grid 
(18x8) 

0 9999 

i 

0 

0. 25802xl0 S + 

0 1560xl0 5 

0.2607xl0 S 

0.99 

i 

0 

0. 35520xl0 S+ 

0. 3430xl0 5 

0 3157x10 s 

0 99 

0 

8 

0.351 i0xl0 5 = 

0 3403xl0 S 

0 3561xlO S 

0 99 

0 

6 

0. 34080xl0 S = 

0 3230xl0 5 

0.3415xl0 5 

0 .99 

0 

4 

— 

0.2831xl0 5 

0.3058xl0 S 

0.99 

0 

2 

-- 

0.2311xl0 S 

0.2449xl0 S 

0.9 

1 

0 

0.7476x10 5+ 

0.7446xl0 S 

0.7483xl0 S 

0.9 

0 

8 

0.5844x10 5= 

0 5718x10 s 

0.5749xl0 S 

0.9 

0. 

6 

0.45360x10 S= 

0.4209xlO S 

0. 4232xl0 S 


+ Happel (1958) 

= Mohan and Raghuraman (1976a) 
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drag coefficients presented in the table*} %, It is thus safe to 
conclude that the 18 x 8 mesh is adequate from accuracy view point 
and computationally economical in creeping flow conditions for 0 3 
< € < 0.9. Indeed, the computational time increases by a factor 
of 6-7 for the 18 x 16 mesh as compared to that of the 18 x 8 
mesh . 

(b> Me?e?d for an Exponential Grid for Dilute Beds: 

A comparative study of a uniform grid and an exponential grid 
has been carried out to demarcate the ranges of their 
applicability. In creeping flow regime, typical results on drag 
coefficient obtained using the aforementioned two grids are 
presented in table 4.9 for Newtonian and power law type fluids and 
for a range of values of « For the sake of comparison, 
analytical results of Happel (1958) for Newtonian fluids and 
approximate results of Mohan and Raghuraman (1976a) for 
non-Newtonian fluids are also included in this table. 
Intuitively, one would expect the two grids to yield comparable 
results for low values of voidage An inspection of Table 4.9 
suggests that the exponential grid consistently yields better 
results than those obtained using the uniform grid under all 
conditions examined here. However, the results from the two grids 
show noticeable differences only for dilute systems (« > 0.9). 
For instance, in the case of € = 0.9, the two results are almost 
identical for both Newtonian and power law liquids, with the 
maximum divergence between the two being only 0.54% in the case of 
n = 0.6. As expected, the performance of the uniform grid 
deteriorates rapidly with the increasing value of voidage. For 
example, at € = 0.9999, the uniform grid provides a very coarse 
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discretization of the domain adjacent to the particle surface and 
gives highly inaccurate results. The use of an exponential grid, 
on the other hand, leads to results which are accurate to within 
1.2% of the corresponding analytical results for Newtonian fluids. 
Similar observations can be made for power law fluids also From 
the aforementioned discussion, it is evident that for e > 0.9, an 
exponential grid system is preferable in order to obtain 
sufficiently accurate results, as compared to a finer uniform mesh 
which will be prohibitively expensive. 

C c) Effect of Physical and Numerical Factors on Rata of 

Convergence of Solution: 

For iterative solutions of highly nonlinear equations (such 
as those of the present study), an important criterion for judging 
the scheme efficiency is the rate of convergence of the solution. 
Usually two criteria, namely, a pointwise convergence for velocity 
and pressure fields, and an overall convergence for the integrated 
quantity such as drag coefficient are applied. Though the 
pointwise convergence is desirable at every node in the flow 
domain, it must be used with a little bit of caution, and no extra 
efforts should be expended in achieving convergence for small 
values of velocity and pressure which are known to undergo wild 
fluctuations during simultaneous solution. Therefore, cut-off 
values for pressure and velocity were specified as 5% of the 
maximum velocity and 7% of the maximum pressure; no attempt was 
made to meet the pointwise convergence criteria below these 
values The results were assumed to have converged if the 
difference in the value of a variable at any nodal point for two 
successive iterations was less than 0.5%. this criterion was 
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sufficient to ensure the convergence upto four significant digits 
in most cases except for a few pressure values Drag force was 
assumed to have converged when no change in the first three digits 
was detected for three to four consecutive iterations. 

Generally speaking, the number of iterations required to 
achieve the pre- determined level of convergence for primitive flow 
variables (u,p.v) and for drag depends upon the complex 
interactions between physical factors such as voidage, Reynolds 
number and power law index and the numerical parameters such as 
initial values, relaxation factor, etc. Consequently, each run is 


di f f er ent 

from 

the other In 

some 

way, 

and indeed it Is 

very 

difficult 

to draw specific 

conclusions 

on this 

aspect . 

Some 

typical 

cases 

are listed 

In 

Table 

4.10 and 

based on 

the 


information given here, the following observations can be made: 

(i) For a given value of Reynolds number and voidage, the 

number of iterations to achieve the desired level of 
pointwise and overall convergence increases with increasing 
non-Newtonian behaviour This finding is in accordance 
with the observations made by Gu and Tanner (1985). 

(ii) For given values of « and n, the number of iterations is 

virtually unaffected by the value of Reynolds number in the 
range employed here; similarly, for fixed values of n and 
Re, the effect of voidage on the rate of convergence is 
insignificant, although there is a weak tendency for the 
number of iterations to increase with decrease in « . 

(iii) In the packed bed range (<s = 0 3 - 0.5) and for highly 

non-Newtonian fluids (n = 0.2 and 0.4), the drag 

coefficient was found to converge between 15 and 20 
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Table A. 10 Influence of various physical and numerical 
parameters on rate of convergence. 
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1 

0.99 

0.8 

0.001 

Trivial 

1.0 

7 

3-4 

2 

0.99 

0.2 

0.001 

Solution for 
e = 0.99, 
n = 0.4, 

Re= 0.001 

1.0/0. 7 

1 1 

7 

3 

0.99 

3.0 

1 

Trivial 

1.0 

5 

3 

4 

0.99 
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Trivial 
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0.99 
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1 

Solution for 
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0.99* 

0.4 

1 

Solution for 
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1.0/0. 7 

8 

5 

8 

0.9 

0.8 

0.001 

Trivial 

1.0 

6 

4 

9 

0.9 

0.8 

1 

Trivial 

1.0 

6 

4 

10 

0.9 

0.2 

1 

Solution for 
e = 0.9, 
n = 0.4, 

Re= 1 

1.0/0. 7 

22 + 

17 

11 

0.9 

0.8 

10 

Trivial 

1.0 

7 

i 

5 

12 

0.9 

0.2 

10 

Solution for 
e = 0.9, 
n = 0.4, 

Re= 10 

1.0/0. 7/0. 5 

nr 

25 

20 

13 

0.5 

0.8 

1 

Trivial 

3.0 

7 

4 

14 

0.5 

0.8 

20 

Trivia 1 

1.0 

4 

3 



15 

0.3 

0.2 

1 

16 

0.3 

0.8 

10 

1 7 

0.3 

0.2 

10 


18 

0.4 

1.0 

10 

19 

C.4 

1.0 

10 

20 

0.9999 

0.8 

10 

21 

0.9999 

0.8 

20 

22 

0.9999 

0.2 

20 


Output 

solution for 
e s 0.3, 
n * 0.4 , 

Re = 1 

1.0/0. 7/0. 5 

Trivial 

1.0 

Output 

solution for 
e = 0.3, 
n = 0.4, 

Re= 10 

3. 0/0. 7/0. 5 

Trivial 

1.0 

Trivial 

0.5 

Trivial 

1.0 

Trivial 

1.0 

Output 

solution for 
t = 0.9999, 
n = 0.4, 

Re= 20 

1.0/0. 7 


+ Results did not converge upto desired level. 
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Iterations whereas the pointwise convergence of (u,p,v) 
required about 20 to 30 iterations However, in some 

cases, even with these many iterations convergence could 
not be achieved In general, the desired drag convergence 
was achieved in about half the number of iterations as 
those required for the pointwise convergence of velocity 
and pressure fields 

(iv) For mild to moderate degrees of non-Newtonian behaviour 
(0.6 £ n :£ 1), the initial guess of a trivial solution for 
pressure and velocity variables led to faster rate of 
convergence than that by using the solution for the case 
having same kinematic (Re) and physical conditions («) but 
a slightly different value of the power law index (n) . 



CHAPTER V 


RESULTS AND DISCUSSION - NEWTONIAN FLOW PAST MULTI-PARTICLE 

ASSEMBLAGE 

In this chapter, results predicted by the free surface cell 
model for steady Newtonian flow past a multi-particle system are 
presented . Pressure , friction and total drag coef f icients , the 
variation of angular and radial velocity components , surface 
pressure and shear stress distributions and streamline plots have 
been obtained f or a wide range of kinematic and physical 
conditions covering 0.001 i Re £ 100, 0.3 < & < 0 9. While the 
above range of voidage represents the complete spectrum of 
mult 1 -part icl e flow situations including fixed, fluidized and 
distended beds , the performance of the free surface cell model in 
the case of approach to a single translating sphere has been 
investigated by using the very high voidage value of <s — 0.9999. 

The numerical results presented herein have been validated by 
means of comparison with the appropriate theoretical and 
experimental results for packed beds, fluidized beds and single 
spheres, available in the literature. The range of applicability 
of the free surface cell model has been identified, and its 
performance has been contrasted with that of its rival, the 
zero-vorticity cell model . Typical streamline patterns are 
presented to elucidate the salient features of the flow field. 



89 


S.l EFFECTS OF VOIDAGE AND REYNOLDS NUMBER ON DRAG COEFFICIENTS 

The complete set of drag results obtained in this study are 
summarized in Table 5.1. The dependence of the total, pressure 
and friction drag coefficients on voidage and Reynolds number is 
shown in Figures 5.1 to 5 3 for representative parametric values. 

For a given voidage, the drag coefficients (pressure, 
friction and total) decrease with Reynolds number. The fact that 
all the three quantities have the slope of -1 on log-log 
coordinates over some Reynold number range indicates that they are 
all proportional to 1/Re in that range. Such behaviour is the 
characteristic of viscosity-dominated flows Indeed, an 

examination of the Happel ’ s analytical solution (discussed in 
Appendix I) reveals that the inverse proportionality between the 
drag coefficients and the Reynolds number is essentially a 
creeping flow feature. Therefore, it appears that for flow 
through multi-particle assemblages the viscous effects remain 
important even in the intermediate Reynolds number range (upto 
about Re ~ 30 - 40). However, beyond a critical value of Reynolds 
number, the drag results begin to veer away from C D Re = constant 
behaviour. The critical value of Reynolds number appears to be a 
function of voidage. For instance, withe = 0.9, this transition 
seems to occur at about Re ~ 10-15, and it takes place at higher 
values of Reynolds number for lower value of voidage. In order 
to understand these trends, it must be recognised that the 
inertial and viscous forces in multi-particle systems are 
functions of voidage. Therefore, the particle Reynolds number 
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Table 5.1: Summary of drag (pressure, friction and total) 
coef f icients for Newtonian fluid flow past an 
assemblage of particles. 


Re 

C DP 

C DF 

C D 

€ = 0.3 




0.001 

0 .4337xl0 7 

0.1210xl0 7 

0.5543xl0 7 

0.1 

0.4333xl0 5 

0 1210xl0 5 

0.5543xl0 5 

1.0 

0.4333xl0 4 

0 1210xl0 4 

0.5543xl0 4 

5.0 

0.8646xl0 S 

0.2418xl0 S 

0.1106xl0 4 

10 . 0 

0.4333xl0 S 

0.1211xl0 S 

0.5544xl0 S 

20.0 

0.2167xl0 S 

0. 6051xl0 2 

0.2769xl0 S 

50.0 

0.8701xl0 2 

0. 2436xl0 2 

0.1114xl0 S 

100.0 

0.4427xl0 2 

0. 1242xl0 2 

0 . 5669xl0 2 


•e = 0.4 

0.001 
0.1 

1.0 

5.0 

10.0 

20.0 

50.0 


0.1438xl0 7 

0.1438xl0 5 

0.1438xl0 4 

0.2876xl0 3 

0.1439xl0 3 

0.7203X10 2 

0.2907xl0 2 

0.1495xlQ 2 


0. 6085xl0 6 
0. 6085xl0 4 
0.6086xl0 3 
0.1217xl0 S 
0.6G91xl0 2 
0.3055X10 2 
0.1245xl0 2 


0.2047xl0 7 
0.2047xl0 5 
0.2047xl0 4 
0 4094xl0 S 
0.2048xl0 S 
0. 1026xl0 S 
0.4152xl0 2 
0.2151X10 2 


100.0 


0.6595x10 
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« = 0.5 


0.001 

0 . 567 3xl(j 

0.1 

0.5673x10* 

1.0 

0.5673xl0 3 

5.0 

0. 1135X10 3 

10.0 

0.5677xl0 2 

20.0 

0. 2845xl0 2 

50.0 

0 .1155xl0 S 

100.0 

0.5992x10 


0.3438xl0 6 

0 . 9111x10** 

0.3438x10* 

0.911 1x10* 

0.3438xl0 S 

0.9111xl0 S 

0.6880xl0 S 

0.1823xl0 S 

0.3455xl0 S 

0.9122xl0 S 

0.1732xl0 S 

0.4577xl0 2 

0.7178x10 

0.1872xl0 2 

0.3887x10 

0.9875x10 


€ = 0.7 


0.001 

0.1144xl0 6 

0.1289xl0 6 

0.2433x10 

0.1 

0.1144x10* 

0.1289x10* 

0.2433x10 

1.0 

0.1144xl0 S 

0.1290xl0 S 

0.2433x10 

5.0 

0.2288xl0 S 

0.2584xl0 S 

0.4872x10 

10.0 

0.1146xl0 S 

0.1299xl0 S 

0.2455x10 

20.0 

0.5768x10 

0.6656x10 

0.1242x10 

50.0 

0.2536x10 

0.2937x10 

0.5293x10 

100.0 

0.1214x10 

0.1665x10 

0.2883x10 


6 

4 

3 

2 

2 

2 



0 001 
0 1 

1.0 

5.0 
10.0 

20.0 

50.0 

100.0 

€ = 0.9999 
0.001 
0.1 

1.0 

5.0 

10.0 

20.0 

50.0 

100.0 


0.2643xl0 5 
0. 2643xl0 3 
0 . 2643xlQ 2 
0.5291x10 
0.2652x10 
0.1337x10 
0.5452 
0.2495 

0.8225xl0 4 

0.8214xl0 2 

0.8273x10 

0.1887x10 

0.1198x10 

0.8491 

0.6196 

0.5123 


0.4803xl0 5 

0.4803xl0 3 

0.4805xl0 2 

0.9709x10 

0.5000x10 

0 2745x10 

0.1504x10 

0.8787 

0. 1784xl0 5 
0. 1787xl0 3 
O.lSOlxlO 2 
0.4261x10 
0.2811x10 
0.2017x10 
0 1338x10 
0.8813 


0. 7446xl0 5 
0.7446xl0 S 
0. 7448X10 2 
O.lSOOxlO 2 
0.7653x10 
0.4081x10 
0.2049x10 
0 1128x10 

0.2607xl0 5 

0.2608xl0 S 

0 . 2628xl0 2 

0.6148x10 

0.4009x10 

0.2866x10 

0.1958x10 

0.1394x10 



ORAG COEF 



0-01 01 1 10 100 
REYNOLDS NUMBER . Re 


Fig. 5-1 Dependence of total drag coefficient on Reynolds 
number and bed voi dage 
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which has been defined without Incorporat ing the effect of 
void&ge, does not represent the true ratio of inertial and viscous 
forces in a dense or moderately dilute assemblage of particles. 
It Is only in the case of extremely dilute assemblages (i.e. 
approaching the limit of a single sphere) that particle Reynolds 
number represents the actual ratio of inertial and vise*, a forces. 
The predominance of viscous force for highly dense assemblages can 
be attributed to narrow pore spacing available for flow. 

At higher Reynolds numbers the versus Re curve becomes 
less steep. This is consistent with the well-known trend for flow 
around single sphere, wherein the slope of vs. Re curve 
decreases in magnitude at higher Reynolds number due to inertial 
effects, until becomes nearly a constant in the turbulent flow 
regime The curious feature, however, is that the decrease m the 
slope of the curve Is largely contributed by the viscous drag 
as shown in Fig. 5.3. 

Figures 5.1 - 5.3 also indicate that for a given Reynolds 
number, the drag coefficients decrease in magnitude as the voldage 
increases. The explanations for this observation are readily 
found if one considers the fact that the available flow area 
reduces when voidage is lowered. In addition to this, the flow 
undergoes expansion and contraction more number of times due to 
the larger solid concentration These factors combined together 
increase the pressure drop and the viscous losses. It is to be 
mentioned here that for steady flow situation, where inertial 
forces are not dominant , the pressure forces are locally balanced 


by the viscous force. 
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An interesting feature, not reported or mentioned in the 

literature hitherto can be clearly seen when the two individual 

components of the drag coefficients are plotted together as shown 

in Figure 5.4 An inspection of this figure suggests that at low 

values of voidage (<e = 0.4), the total drag coefficient draws a 

major contribution from pressure drag. However, as the voidage 

increases, the contribution of friction drag gradually increases 

and becomes nearly equal to the pressure drag at « = 0.65 or so, 

and finally it dominates over the pressure drag (« = 0.9) The 

above discussed trend implies that the contribution of pressure 

and viscous drag to the total drag is a function of bed voidage. 

An order of magnitude analysis (see Appendix- 1 1 ) shows that 

3 

pressure and viscous forces vary respectively as 0(1/6 ) and 
2 

0(1/6 ) , where 6 is the radial gap between the particle and the 
cell surface. It is expected that the pressure and friction drag 
coefficients, which are integrated quantities arising from the 
surface stresses should show atleast qualitatively a similar 
dependence on 6. It can be proved from such an order of magnitude 
estimates that 



For 



FRICTION AND PRESSURE DRAG COEFFICIENTS 


Fig. 5 



4 Effect of Reynolds number end bed voidage on 
relative magnitudes of pressure and friction drag 
coefficients . 
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Nj shi o.ura a'td lahii (1980a) asserted that the pressure drag 
coeft icient , C^p is independent of the value of Reynolds number 
for Re > 10 for all values of voidage m the range 0 1 :£ e £ 0.3. 

The observation is not in conformity with the present numerical 
r esul ts . 

In another study Leclair and Hamlelec (1968) found that the 
C 

* . DP . 

ratio, ts increases with increasing Reynolds number, after the 
DF 

initial range in which it remains constant The present numerical 

results do not display that trend For a given value of bed 

C 

. DP 

voidcge, ts — remains constant initially upto a certain value of Re 
DF 

and then starts decreasing with Re. It must be borne in mind, 
that the earlier authors have not explained their results 
physically, whereas, the present results are found to be 
physically consistent with the order of magnitude estimates. 


5.2 SURFACE PRESSURE AMD SHEAR STRESS DISTRIBUTION 

It is Intuitively expected that the study of pressure and 
surface shear stress distribution over the sphere surface will 
facilitate further rationalization of some of the phenomena 
discussed In previous section. The surface values of both 
pressure and shear stress have been re-normalized with respect to 
a characteristic viscous stress (rj U /R D ) and are shown m Figures 
5.5 to 5.16, for various combinations of parametric values to 
elucidate the effects of voidage (e) and Reynolds number. From 
the results shown in these figures, the following observations can 
be made: 

(i) At creeping flow (Re = 0.001), the numerically predicted 


Dimensionless surface shear stress 



Fig. 5.5 Effect of Reynolds number, on the variation 
of nondimensional surface shear stress T g /(>? 
with angular position, for « = 0.4. 


Dimensionless surface shear stress 


1 



Angular position ,8 


Fig. 5.6 Effect of Reynolds number, on the variation of 
nondimensional surface shear stress t s /(t) U^/Rp) 
with angular position, for e = 


0.7. 




Dimensionless surface shec. stress 
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0 30 60 90 120 120 180 

Angular position ,0 

Fie 5.7 Effect of Reynolds number, on the variation of 

nondimensional surface shear stress r / (rt U /R„) 

s eo P' 

with angular position, for « = 0.9. 
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surface shear stress varies sinusoidally with angle at all 

voidages, in conformity with the Happel’s analytical 

solution Moreover , the shear stress exhibits only 

marginal deviation from the creeping flow results even upto 

the Reynolds number of 20 at a low voidage (« = 0.4). This 

is a consequence of large contact area existing between the 

solid and the fluid phases when « is small At higher 

voidages (« = 0.7 and e = 0.9) however, there is noticeable 

divergence from the creeping flow results even at Re =10. 

In the intermediate Reynolds number regime, the surface 

shear stress distribution is asymmetric about the mid plane 

(0 = n/ 2) with increasing skevedness as the bed voidage 

increases. The maximum shear stress shifts towards the 

rear region (90° < & < 180°) for such a situation. This is 

believed to be resulting from fluid Inertia, which causes 

the fluid to move with nearly the particle-velocity in 

front of it. It is evident from Figures (5.5 - 5.7) that 

in the range 0.001 < Re < 20 and Q . 4 < <= < 0 9, the voidage 

has a stronger influence on the dimensionless shear stress 

than the Reynolds number. The weak effect of Reynolds 

number is explained in terms of the fact that for the 

multiparticle flow situation, the role of inertial force is 

minor even upto Re = 20. As regards voidage, the claim 

2 

that shear stress varies approximately as 0(l/«5 ) is 
supported by trends shown in Figures . 

(ii) The pressure profiles (Figures 5.8 - 5 10) vary according 
to the cosine of the angle at low Reynolds number (Re < 1) 



Non dimensional surface pressure 


1 



Void troction (£)*0 4 


Angular position, 6 


Fig. 5.8 Effect of Reynolds number, on the variation 
nondimensional surface pressure P 8 /(n U «/ R p) 
angular position, for « = 0 4. 


Nondimensionol surface pressure 





5.9 


Effect of Reynolds number, on the 
nondlmen.ionsl surface pressure p,/C» 
angular position, for & 


variation of 
U^/Rp) with 



Nondimensional surface pressure 



Effect of Reynolds number, on the variation of 
surface pressure p g / (?) with angular 
position, for c = 0.9. 


Fig. 5.10 
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for all voidage® considered here. This is In accordance 
vith the analytical results of Rappel . However, unlike the 
surface shear stress profile, the surface pressure at 
i nt ermediat e Reynolds number range exhibits larger 
deviation from creeping flow results even for low voidages . 
Thus, it appear® that the inertial effects more directly 
Influence the pressure than the shear stress At higher 
voidages, the veering away of the surface pressure 
variation from corresponding creeping flow pressure profile 
Is even more obvious This is hardly surprising because 
pressure is more sensitive to changes in voidage thar the 
shear stress, in accordance with the order of magnitude 
estimates discussed earlier in connection with the pressure 
and friction drag coefficients. Indeed, for Re = 20, when 
■the voidage is changed from £ = 0.4 to 0.9, the maximum 
shear stress reduces only by a factor of 7 or 8, while the 
stagnation pressure reduces by a factor of 50. 

(iii) The relative magnitudes of the surface pressure and the 
surface shear stress at three representative voidages (« = 

0. 4, 0.7 and 0 9) and two Reynolds numbers (Re = 0.001 and 
20) are shown in Figures 5.11 to 5.16. It is observed that 
at low voidage (e = 0.4), the dimensionless pressure 
variation over almost the whole particle surface is much 
larger in magnitude than that of the dimensionless shear, 
save for a small part of the surface in the rear of the 
sphere. As the voidage increases, the trend reverses, 

1. e., the front half of the sphere is dominated by pressure 



Nondimensional Surface pressure 8. shear stress 





Nondimensional surface pressure 8. shear stress 


Surface pressure 

Shear stress 


30 60 90 120 

Angular position, 6 


S.12 Relative magnitudes of nondxmensiona, surface Rressur, 

P /(„ U /R p ) and ahear etreas rjln VV’ for ‘ ' 
s <40 r 

and Re = 0.001. 


Non dimensional surface pressure & shear stress 



Non dimensional Surface Pressure And Sheer Stress 



ANGULAR POSITION , 6 


Fig - 5.14 Relative magnitudes of nondimensional surface 

pressure p /(tj U /R D ) and shear stress r /(rj U /R e ), 

8 Jr s ® r 

for « = 0.4 and Re = 20. 


Non dimensional Surface Pressure And Shear Stress 



ANGULAR POSITION , 6 

Fig 5 Relative magnitudes of nondimensional surface 

pressure p # /(r> U^/R p ) and shear stress t ,/(» VV 
for « = 0.4 and Re = 20. 


Nondimensional Surface Pressure And Shear Stress 



ANGULAR POSITION , 0 

Fig . 5.15 Relative magnitudes of nondimensional surface 

treasure p e /(r> U^/Rp) and shear stress t g / ( t, U^/Rp), 
for « = 0.7 and Re = 20. 




Fig. 5.16 Relative magnitudes of nondimeneional surface 

pressure Pf> /( n U^/Rp) and shear stress t J (» U^/Rp). 
for « = 0.9 and Re = 20 
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and the rear portion by surface shear stress. This becomes 

increasingly evident as the bed voidage increases. 

Finally, with the increasing dilution of the assemblage, 

the ratio (P /t ) decreases from about 9 to 2 for the 
max max 

creeping flow range (Re = 0 001) At « = 0.9, the 

variations in the surface pressure and shear stress becomes 
comparable so that the relative drag contributions are also 
of the same order 

The effect of change in Reynolds number on the relative 

magnitudes of the local values of surface pressure and 

shear stress is shown in Figures 5.14 - 5 16 wherein the 

results are plotted for Re = 20 The trends are 

qualitatively similar to the creeping flow range In the 

packed bed range (s = 0 «), the surface pressure dominates 

over the corresponding surface shear stress over a major 

part of the sphere surface (0° £ & £ 110 ). However, t° r 

distended beds (c = 0.9), the pressure dominated region is 

reduced (0° i£ 9 * 60°) while the > 100°) region is 

clearly dominated by surface shear, with a transition zone 

in between these two parts For Re = 20, the value of 

P /t falls from 11 to 2, a drop comparable to that m 

max max 

the creeping flow region 

5.3 ANGULAR AND RADIAL VELOCITY PROFILES 

It is generally realized that the local values of surface 
shear and pressure are directly determined by the local v 
gradients prevailing in the flow field, in particular, gr 
of the angular velocity component 


F igurts 5 17 and 5 IS show the typical 


variation 
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of 

normal iced for e = 0 4 and 0 9 in the creeping f ov regime 
(Re = 0 001 and 1) An examination of the results plotted in 
these two figures suggests that complete fore and aft symmetry 
prevails as the values of corresponding to & — 30° and 150°, 
and & = 60 and 120° coincide with each other. The results for 
the remaining values of voidage (not shown here) also conform to 
this behaviour. At such low values of Reynolds number, the 
angular velocity seems to change its direction close to the sphere 
surface, at r ~ 1.02 for « = 0.4 and at r ~ 1.033 for « = 0.9. 
Furthermore, regardless of the value of Reynolds number, it would 
be reasonable to expect steeper velocity gradients in low voidage 
multi-particle systems than those for dilute systems (See 
Appendix-I) and this assertion is indeed corroborated by the 
results shown in Figures 5.17 to 5 20. 

The effect of Reynolds number on the angular velocity 
profiles is brought out in figures 5.19 and 5.20. Clearly, the 
flow is no longer symmetrical (also revealed by surface shear 
stress profile), and more significantly, the change of direction 
of the angular velocity now takes place at varying radial 
distances depending upon the angular coordinate (#). The angular 
velocity gradient being maximum in the vicinity (some where on the 
right) of & = 90° plane, the shear stress is expected to go 
through a maximum value at the same location on the sphere 
surface. This is approximately borne out by shear stress profiles 


also . 



Dimensionless angular velocity 


116 



Fig . 5.17 


Variation of nondlmens i onal angular 
radial direction for <e = 0.4; Re = 


velocity in 
0.001 and 1.0 


Dimensionless angular velocity 


Fig 



5.18 Variation of nondintensional angular velocity in 
radial direction for € = 0.9; Re = 0.001 and 1.0 





ro=3 


118 



Fig. 5.19 Variation of nondiutens ional angular velocity 
radial direction for <= = 04 and Re = 20 








Dimensionless radial velocity 


0'6r 


€ = 0-9 
R € =0-001 

© 

h 0 

o 30 
o 60 
• 120 
a 150 

A 180 


- 0*6 


-1-0 J 


Dimensionless 
radial distance- 


«•- 5 - 21 in 





Dimensionless radial velocity 


21 



Fie 


S 22 


Variation of nondineneional 
radial direction for « = 0-* 


radial velocity 
and Re = 20. 
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Though It is generally recognized that in the type of flow 
being considered here, the angular velocity plays a more important 
role than the radial velocity, it is together informative to 
provide some representative results on the radial velocity 
profiles. Owing to the Imposed boundary conditions of no slip and 
the free stream velocity at the two boundaries, the normalized 
radial velocity varies between 0° and ± cos O m the radial 
direction Typical variation of with r for a range of values 
of & and for € = 0.9, is shown in figures 5 21 and 5.22 for two 
different values of Reynolds number. The results for other values 
of voidage also show qualitatively similar trends. 

5.4 STREAMLINE PATTERNS 

Streamline plots provide a qualitative picture of the flow 
patterns occurring in the flow field, and also facilitate the 
identi f ication of points of separation, flow reversal, etc. 
Figure 5.23 and 5.24 (drawn to the scale) depict the streamlines 
for two values of Reynolds number (0.001 and 20) for « = 0.9. 
Evidently, complete symmetry exists at Re = 0.001 (Figure 5.23) 
while the analogous plot for Re = 20 (Figure 5.24) shows that 
asymmetry has set In. Streamline patterns for other values of 
voidage also show similar trends except that lower the value of e , 
higher is the value of Reynolds number upto which the symmetry 
persists. However, no separation in flow was observed under any 
combination of conditions of Re and c studied herein. 

After discussing the detailed flow results obtained here, It 
Is desirable and instructive to examine their utility through 
comparisons with the available experimental results for dense and 



^ = 0-654 



0-6587 
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dilute assemblages Such comparisons facilitate to identify the 
range of applicability of the free surface cell model. 

5.5 COMPARISON OF PREDICTED RESULTS WITH EXPERIMENTAL DATA FOR 

PACKED BEDS 

The results presented in the previous chapter indicated that 
excellent correspondence exist between the present results and 
analytical results of Happel (1958) for a creeping flow situation. 
Happel (1958) in his paper made the comparisons of the 
free-surface cell theory results in the Stokes flow range with the 
results of other theoretical analyses and the experimental data 
then available in the literature. The close agreement of theory 
with the well known Carman-Kozeny or Fair-Hatch equations which 
are widely used to correlate experimental data on packed beds, was 
emphasised. The average deviation of Happel ’ s results from the 
prediction of Carman-Kozeny relation (with Kozeny constant, k = 
5.0) over the range of voidage € = 0.2 to 0.7, was reported to be 
8%. The maximum deviation which occurred for the voidage, « = 

0.7, was found to be 13%. Therefore, the validity of the free 
surface cell model for the creeping flow of a Newtonian fluid 
through a packed bed is fairly well-established The objective of 
the present investigation is to identify the applicability of the 
free surface cell model In the intermediate Reynolds number 
regime . 

In the Intermediate Reynolds number regime, the Ergun 
equation (1952) purports to represent most of the available 
experimental pressure drop data on flow through packed beds with 

The Ergun equation is written as. 


reasonable accuracy. 



FRICTION FACTOR . 


Fig 



.25 Comparison of present ™^lcal ^^ ul ^ ow Vi ^ roueh 
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f bk = 150 ~Ee ) + 1 75 (5 1) 

As mentioned earlier in Chapter IV, El-Kaissy and Homsy (1973) 
extended the analysis of Happel (1958) to the intermediate 
Reynolds number range by carrying out a regular perturbation 
analysis. A comparison between their and present results in terms 
of friction factor was shown in Table 4.6. In the same table the 
values of f^^ calculated using the Ergun’s equation are also 
included. The comparison of the present results with those of 
Ergun is shown in Figures 5 25. Evidently, the present results as 
well as those of El-Kaissy and Homsy (1973), are in reasonable 
agreement with the predictions of Ergun’s equation upto about Re ~ 
20 In the range 0.3 S « i 0.5 Beyond this value of voidage, the 
correspondence between the two deteriorates rapidly and the 
divergence becomes almost 100% at Re = 100. If it is assumed here 
that the Ergun equation faithfully represents the experimental 
data, it is then obvious that the free surface cell model, 
originally developed for creeping flow conditions, breaks down at 
about Re ~ 20. This is possibly due to appreciable convective 
effects and it is no longer possible to retain the spherical 
symmetry of the cell. In view of this, the utility of the results 
presented by Fukuchi and Ishii (1983) and Ni shimura and Ishii 
(1980a) for the range 10 < Re < 10 5 is questionable. 

5.6 COMPARISON OF PREDICTED RESULTS WITH EXPERIMENTAL DATA FOR 

FLUIDIZED BEDS: 

In the theoretical modelling of flow through fluidized beds, 
considerable efforts have been spent on the prediction of minimum 
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fluidization velocity and establishing the bed-expansion 
characteristics . Based on experimental evidence, some researchers 
developed equations relating the relative velocity ratio of the 
assemblage to the bed void fraction. Here the relative velocity 
ratio refers to the ratio of the actual velocity of the assemblage 
U^, to the terminal free fall velocity of the single sphere . 
The pioneering work in this direction concerning a homogeneous 
fluidized bed was by Richardson and Zaki (1954). Several other 
investigators such as Lewis and Bowerman (1952), Letan (1974), 
Garside and Al-Dibouni (1977) followed the same approach and 
proposed somewhat different forms of equations. These 
correlations have been summarised by Couderc (1985), Khan and 
Richardson (1985), etc. 

The drag results of the present study were used to provide a 
general fluidization correlation which is more useful for design 
and analysis (Zenz & Othraer , 1960) by replotting the drag 
coefficient variation with Reynolds number and porosity data in 
terms of (Cp Re 2 ) 1 ^ 3 against (Re/Cp) 1 ^ 3 with porosity as a 
parameter . Such a plot is presented in Figure 5-26 The plot is 
useful as it essentially relates the particle velocity particle 
diameter with physical properties as fixed parameters, i.e. 


(C D Re 2 ) 1/3 = 


4 ^ C V P) 6 ,1/3 
P ' 3 


d { ’ § } 


r> 


i/o 3p 

(Re/C D ) ' = < 4y) (p>T 5 


1/3 


(5.2) 

(5.3) 


The relationship between the relative velocity ratio, (U^/l^) and 
the porosity («) can be obtained simply by arawing a line parallel 
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to (Re/Cp) ^ axis through the particular particle diameter or 
free settling Reynolds number desired The variation of relative 
velocity ratio over a range of voidage (0 4 to 1.0) at several 
di f f event Reynolds numbers ( Re^ = 0 1 , 10, 20, 50, 100) are shown 
in Figures 5.27 to 5.31. On these figures, the predictions of 
correlation by Richardson and Zakl (1954a) which is believed to 
repr esent r easonably most of the experimental data are also shown 
for compar Ison purpose. In order to provide more comparison with 
experimental data, the predictions of the popular correlation by 
Gars Ide and Al-Dibouni (1977) are also shown in these figures 

For low Reynolds number (Re^ < 0.1), the present numerical 
results lack good agreement with the predictions of Richardson and 
Zaki as well as Garside and Al-Dibouni . At the lowest end of 
fluidized bed voidage, « = 0 4, the present results are found to 
deviate from those of Richardson and Zaki by 17% which Increases 
to about 30% at € = 0 5. At still higher voidage (« > 0.5) the 
deviation Increases to 50% or more. Earlier Happel (1958) made 
exhaustive comparison of free sut face cell model results for 
creeping flow with large number of experimental results of several 
inves t Igat or s for fluidized beds. He also observed the lack of 
agreement between theoretical predictions and experimental 
results. One of the reasons for such disagreement could be the 
agglomeration and circulation effects which will result in 
increased value for relative velocity. It is to be mentioned here 
that the voidages measured in experiments are average quantities 

f 

and the actual distribution of particles may not be uniform, 
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whereas, the free surface cell model assumes a uniform particle 
distribution. 

At intermediate Reynolds number (1 i Re i 100) the agreement 

o 

between the present numerical results and predictions of 

Richardson and Zaki (1954) considerably improves. At Re Q = 20, the 

two results show excellent agreement and except at « = 0 4 where 

the deviation between the two results is 25%, at all other 

voidages the two results differ by less than 6% At Re = 50 the 

o ’ 

deviation between present results and those of Richardson and Zaki 

again increases At low voidages (« = 0 4 and 0 5), it is 

observed to be about 50% which is found to be limited to 20% at 

higher voidages (0.9 < <£ < 0 6 ) At Re = 100, the deviation 

o 

between the two results are found to be more than 60% at lower 
voidage and around 40% for higher voidage. There are several 

possibilities for the observed trend of results; (a) The free 
surface cell model breaks down beyond Re^ > 20, (b) Secondary 

effects take place in actual experiments 

5. 7 COMPARISON OF PREDICTED RESULTS WITH EXPERIMENTAL RESULTS FOR 
SINGLE SPHERES* 

The behaviour of free surface cell model in dealing with the 
extremely dilute assemblages has been identified by comparing the 
present results with the prediction of the empirical correlations 
to estimate the resistance exper i enc ed by a single sphere in its 
motion through a fluid. It is well known that an analytical 
solution has been so far possible only for the creeping flow past 
a sphere (Re < 0.1). For intermediate Reynolds number, numerical 
solutions (Hamielec et.al., 1967; RImon and Cheng, 1969) are 
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available but they are limited to Re £ 400 Numerous empirical 

correlations based on experimental data have been proposed by 
investigators to cover a wider range of Reynolds number. The 
effect of boundary has also been extensively studied These 

empirical correlations have been summarised by Clift et.al. (1978) 
in their excellent treatise. More recently Khan and Richardson 
(1987) reviewed various empirical correlations available in 

literature. One of the correlations given by Lapple and Shephard 

has also been reproduced in Chemical Engineers Handbook (Perry and 
Chilton, 1973). 

In a theoretical modelling, one would intuitively expect that 
when the cell radius becomes infinitely large the cell model 
should approach the limiting case of a single sphere moving in an 

Infinite expanse of fluid. However, it Is virtually impossible to 

simulate these conditions in numerical work owing to the computing 
requirement becoming prohibitively exorbitant. Thus, one has to 
arbitrarily fix the cell boundary at a sufficiently large distance 
to simulate the conditions corresponding to that of a single 
sphere. In the present study, a very high value of voidage O = 
0.9999) has been used to investigate the limiting behaviour of the 
free surface cell model. The above value of <£ refer to the cell 
to particle radius ratio of 21 47. 

From the available evidence on wall effects ( Cl 1 ft et.al. 
1978) for spher es-in-tube configuration, it is known that the wall 
effect decreases with the Increasing Reynolds number One expects 
at least qualitatively a similar effect of cell radius in case of 
sphere- in- spher e configuration of the spherical cell model. 



DRAG COEFFICIENT 



REYNOLDS NUMBER , Re 


Fig. 5.32 


Comparison of present Ji® curve 
(€ = 0.9999) with Lapple Shephard curve 
« 1 . 
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The present numerical results (for « = 0 9999) are compared 
In Figure 5.31 with the widely used Lappl e-Shephard curve in the 
range 0.001 i Re i 100. The results are in good agreement (within 
5% or so) upto about Re i* 30 but begin to diverge appreciably at 
higher Reynolds number values It is, however not at all clear 
whether such a poor correspondence (for Re > 30) is due to 
inaccuracy of numerical results or due to degree of uncertainty 
inherent in Lappl e-Shephard equation or due to the combination of 
both. Uith regards to the creeping flow regime, it was mentioned 
m Chapter IV that present drag results for <s = 0.9999 are m 
close agreement (within 1.3%) with the analytical results of 
Happel (1958). It is also to be mentioned here that Gu and Tanner 
(1985) studied the creeping motion of a sphere in the power law 
type non-Newtonian fluids Though they have not stated it 
explicitly, their flow configuration and boundary conditions are 
equivalent to the one used in the present study. They have given 
a formulae which allows the estimation of the extent of wall 
effects (for sphere-in-sphere configuration) on drag coefficients 
for the creeping motion of a sphere in Newtonian fluid. This is 
expressed as 


'D 


= c Dst [1 ~ ^ it + 0 ( ir )5 ] 1 


(5.4) 


are the 


where C Dgt is the Stokes drag ( c Dat “ 24/Re) R c and R p 
cell radius and particle radius respectively. For € = 0.9999 and 
Re = 0.001, the present result agrees with that of Gu and Tanner 
within 1.3%. In fact, the drag coefficient evaluated from above 
expression is exactly equal to the drag coefficient obtained from 
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the analytical expression by Happel (1958). In the light of 
foregoing discussion, it Is evident that free surface cell model 
satisfactorily predicts the momentum transfer In a dilute 
assemblage of particles upto Re ^ 20 and this may safely be taken 
as the upper limit for applicability of free surface cell model. 

Before concluding this chapter, it is felt desirable to 
compare the performance of present model with other spherical 
model, namely, the zero-vort ici ty cell model in the following 
section. 

5. 8 COMP AM; SOM WITH ZERO VORTXCITY CELL MODEL 

It Is of interest to assess the relative performance of the 
free surface and zero vorticity cell models over a range of 
Reynolds numbers. Numerical results for the latter model have 
been presented by Leclair and Hamielec (1968). The predictions of 
the two models for packed beds (0.3 ^ « £ 0.6) have been compared 
In Table 5.2. The predictions of Ergun equation are also included 
in this table. In the packed bed range there are only two values 
of voidage which are common to both studies. An examination of 
Table 5.2 shows that the free surface cell model performs 
consistently better than the zero vorticity model in the creeping 
flow region (Re < 10) whereas, the latter is a better predictor in 
the range Re > 10 or so. Furthermore, it is clearly seen that the 
zero vorticity cell model always yields the values of friction 
factor which are higher than those based on free surface cell 
model . This is mainly due to the exchange of energy with other 
cells In the case of zero vorticity model. This conclusion Is in 



Table 5.2: Relative performance of free surface and zero vorticity cell 
models against the predictions of Ergun's equation. 
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line with the findings of others ( e g. Happel and Brenner, 1965; 
El-Kaiasy and Homsy , 1973) 

The performance of the two models in predicting the relative 
velocity ratio (Superficial velocity/Free fall velocity of a 
single sphere) for flow through fluidized beds at low Reynolds 
number has been compared in Fig. 5.28 The results of free 

surface cell model are found to be in close agreement with the 
predictions of Richardson and Zaki correlation. 

5.0 CONCLUSIONS 

Detailed comparison with the analytical solution of Happel 
suggests that the present numerical results are accurate to within 
1.51: in the creeping flow regime It also suggests that the free 
surface cell model, owing to strong convective effects, does not 
provide a good description of flow through multi-particle systems 
for Re £ 30. The analysis of detailed results has also clearly 
demonstrated that as far as the drag values are concerned, the 
approximation of creeping flow (CpRe = constant) holds good upto 
about Re is 10-20 depending upon the value of bed voidage. This 
conclusion is also corroborated by the surface pressure and shear 
stress variations in the same range of Reynolds number values. 
Streamline plots lend further support to this conclusion. The 

pressure field deviates most from the creeping flow results among 
all the flow variables as Reynolds number is increased, although 

this does not affect the drag behaviour significantly. The total 

• „ ..A^rtlons of contributions from 

drag coefficient draws varying prop 

, . respectively depending upon the 

pressure and surface shear stress resp 

In low voidage ayatena, the drag ia do.inated 


value of voidage. 
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by pressure but the latter contribution decreases as the voidage 
increases This phenomenon has been explained qualitatively via 
the angular velocity profiles In the limiting case of a single 
sphere, notwithstanding the complexity of wall effects, the 
present results in the creeping flow region show an excellent 
agreement with the existing analytical solution, however, the 
correspondence is only moderate for Re > 20 

The relative velocity ratio data for fluidized beds show lack 
of agreement with the predictions of Richardson and Zaki 
correlation in the creeping flow regime However, the agreement 
between the two improves with incer easing Reynolds number, at Re = 
20, the deviation between the two results for all voidages 
corresponding to fluidized bed system (0 8 > -e > 0 4), is found to 
be less than 8%. Uith further increase in Reynolds number, the 
two results begin to diverge again 

Based on the comparison with the widely used Ergun equation 
for flow through packed beds, it appears that the free surface 
cell model is a better predictor upto about Re < 10 and beyond 
this value, the zero vorticity cell model is superior, thereby 
making them complimentary to each other. 



CHAPTER VI 


RESULTS AND DISCUSSION - NON-NEWTONIAN FLOW PAST MULTI-PARTICLE 

ASSEMBLAGES 


Steady flow of an incompressible power law type non-Newtonian 
fluid through an assemblage of rigid spherical particles has been 
invesigat ed . Comprehensive results on the detailed structure of 
flow field as well as the macroscopic parameters such as drag 
coefficients have been obtained over a wide range of fluid 
behaviour ( 1 < n < 0.2), bed voidage (0.3 < « < 0.9999) and flow 
conditions (20 £ Re £ 0.001). The variation of the total and the 
individual drag coefficients (C^, and C^p) under the combined 
influence of the non-Newtonian charact eristics of the fluid and 
the concentration of assemblage are highlighted with the help of 
the corresponding stress and velocity profiles existing around a 
typical particle. Although, drag coefficients and flow field 
results have been computed numerically at several values of 
Reynolds number in the range 1 ^ Re ^ 20, only the drag behavior 
at Re = 20 has been discussed with the aid of the surface stresses 
in section 6 . 2* for the sake of brevity. Finally* the numerical 
results reported herein are compared with the appropriate 
theoretical and experimental results available in literature for 
packed beds, fluidized beds and flow around single spheres. 

The objective of the present investigation has been to 
elucidate the effects of n and * on the flow field and dra 8 
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coefficients for both creeping and intermediate Reynolds number 
flows. The flow field at intermediate Reynolds number is 
Inherently different from the creeping flow field due to the 
presence of inertial forces. Therefore, in order to distinguish 
the effects of inertial forces, the results for low and 
intermediate Reynolds number flows have been discussed here in two 
separate sections . 

0.1 LOW REYNOLDS NUMBER FLOW 

In this section, the results are presented for only one value 
of Reynolds number, i e . Re = 0.001; this value of Re is 
considered to be sufficiently small so as to render the inertial 
effects negligible. As stated earlier, the objective has been to 
Identify the combined influence of flow behaviour index and bed 
voidage on flow field and drag coefficients. 

0.1.1 Drag Phen&onte&raa 

Friction , pressure and total drag coefficients for the range 
of voidages («e = 0.3 to 0.9999) and power law indices (n = 0.2 to 
0.8) are plotted in figures 6.1, 6.2 and 6.3. For reference 
purpose the drag coefficients are also summarized in Table 6.1. 

Uith reduction in n, the viscous resistance weakens owing to 
the shear thinning behaviour of the fluid and this is manifested by 
a decrease in friction drag coefficient. The variation of log C^p 
with n is almost linear from n = 1 to 0.5 and appears to obey a 
relationship of the form C^p oc x , upto about n m 0.5. However , 
as the value of n drops below 0.5, the numerical results veer away 
from the postulated power law behaviour. By analogy with the 
power law fluid model, one would intuitively expect the quantity x 
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Table 6.1 Drag Coefficients for Creeping (Re = 0.001) motion 
of an assemblage of spheres in power law fluid 
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Friction drag coefficient , Cdf 


Fig. 6 



1 Variation of friction drag coefficient with power 
law index and voidage , for Re = 0.001 



Pressure drag coefficient , C 



Fig. 6.2 Variation of pressure drag coefficient with power 
law index and voidage, for Re = 0.001. 


Total drag coefficient , Cp 
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Fig- 



Variation o£ total drae co.fficiant with povar lav 
voidaaa, for Ka - 0-OQ1. 
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to be some sort of average shear rate on the particle surface. An 

order of magnitude analysis (outlined in Appendix — II) indeed 

suggests that the value of x is approximately proportional to the 

average surface shear rate of a corresponding Newtonian fluid flow 

situation with the same values of Re and « 

A further examination of figure 6 1 reveals a decrease in the 

value of C^p with increasing voidage (i e inter-particle 

separation) for a given value of n. It can be readily shown that 

in the case of dense assemblages, the average shear rate on the 

2 

particle surface varies as 0(1/6 ), for a fixed value of n. Here 

6 refers to the radial gap between the particle and cell surface. 

It is also observed that the slope of log vs. n curve 

decreases with voidage. This lends further credence to the power 

law type of relation suggested above. 

In the case of an extremely dilute assemblage (« i* 0.9999 ; 

R c /Rp ~ 20), the dependence of C^p is characterised by an initial 

increase (0.6 < n i 1) followed by a downward trend with 

increasing degree of shearthinning behaviour (n ■& 0.6). It is 

presumably so due to the fact that C^p is predominantly determined 

by the shear stress prevailing m the vicinity of & - nf 2 plane. 

The observed trend of C-- with respect to n as e ■» 1.0 is a direct 

DF 

consequence of shear stress behaviour around 6 - n/Z and this 

aspect is discussed in detail subsequently. 

Now turning our attention to the behaviour of the pressure 
drag coefficient (C Dp ) shown in Figure 6.2, it Is seen that C Dp 
varies with n and <e In a qual i tat ively similar manner to that of 
Cpp. This is hardly surprising, since in creeping flow regime, 
the pressure and viscous forces locally balance each other 
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everywhere in the flow field In the case of dense assemblages 
(low values of €) , the pressure drag coefficient contributes more 
towards the total drag coefficient, this trend continues all the 
way upto € = 0 99. However, as the limit of single sphere is 
approach^* this trend reverses, i e , the total drag coefficient 
draws a larger contribution from the frictional drag than from the 
pressure drag This phenomenon can be clearly seen in Figure 6 4 
where the ratio (C^p/C^p) r s plotted as a function of n and « . 
The switch over in the behaviour of C^p and C^p with variation in 
n in the limit of flow over a single sphere is possibly the reason 
for the drag coefficient attaining a maximum value around n ~ 0.6 
or so, as has been reported by Crochet et.al., (1984) and Gu and 
Tanner (1985) . 

The fact that the total drag coefficient draws a larger 
contribution from the pressure drag in the case of dense 
assemblages, can be explained by the following simplified 
analysis. The momentum equation on the particle surface in & 
direction can be written as 


1 d , . _ 1 dp 

r dr < ‘ rT rd ) “ r S& 


( 6 . 1 ) 


Using the boundary conditions (3.20a, b) of the present 
formulation, the surface shear stress can be simlified as: 


T 


r=l 






dr ' 1 r=l 


( 6 . 2 ) 


Now, the orders of magnitudes of and p appearing in equation 

(6.1) can be shown (See Appendix-II) to depend upon <5 as: 


T 


r# 


pL 


r = l 


oc 


« l/<5 2n 

l/6 2n+1 


(6.3) 

( 6 . 4 ) 
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to be applicable in the range 0 . 5 i n £ 1.0 for the total drag 
coefficient also. This feature could be helpful in testing the 
rheological properties of a fluid vis-a-vis the power law model 

The total drag behaviour in the extremely dilute range 
(« i 0.9999) which is close to the case of a single sphere, is 
shown in Figure 6.5. It is to be noted that at this voidage, the 
pressure and friction drag behaviour with respect to n are not 
similar. For instance, the pressure drag monotonically increases 
with decreae in n while the friction drag shows an initial mild 
increase and then a drop in the value beyond n at 0.5. The 
combined effect of such variations of C ntJ and C n „ leads to an 
initial rise in the value of with its maximum occurring around 
n a* 0.6 followed by a sharp reduction for n < 0.6. It is to be 
mentioned here that for a single sphere situation the initial rise 
of total drag is predominantly contributed by the pressure drag 
component while the subsequent drop in its value is solely due to 
friction drag component. In Fig. 6.5, the results of several 
other Investigators available in literature are also plotted with 
the present results for the sake of comparison. A moot point to 
be noted is the wide deviation existing between the predictions of 
different investigators. This has already been illustrated 
earlier by Chhabra et.al. (1980). The present numerical results 
are closest to the lower bound predictions of Cho and Hartnett 
(1983). 

6.1. 2 Surface; Stresses 

It is generally acknowledged that the values of the 
integrated quantities such as drag coefficients, are determined by 
the stresf fields prevailing around the pat tides, although exact 
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Fig 6.6a Variation of non-diineneional surface shear stress 
(t /pM 2 ) with angular position for various va ues 

of power law index and « = 0-«- 


Non-dimensional surface shear stress, r c x 10 
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Variation of surface shear stress ( T ^ 
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fr = 0.4 and 0 2) and = 0.4. 
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correspond enc e may not be exhibited between the surface stresses 
and drag. One can thus possibly attempt to explain the dependence 

oi ^DF' ^DP an< * on n anc * € examining the pressure and shear 
stress distribution over the particle surface. Though extensive 
results embracing wide ranges of values of n and «e have been 
obtained In this study, only representative stress profiles are 
discussed in greater detail here for the sake of brevity. 

Figures 6.6a and 6.6b display the variation of surface shear 
stress with the power law index (n) for a constant value of 
€=04. As expected, the shear stress exhibits complete fore and 
aft symmetry for all values of n studied here, and also it varies 
in a sinusoidal manner for the Newtonian case. However, as the 
degree of shearthinning gradually increases, the surface shear 
stress shows an Increasing tendency to flatten out, though still 
being symmetrical about & = 90° plane. This is solely 
attributable to the combined effects of the variation of local 
viscosity and shear rate. The local viscosity decreases from the 
axis upto the mid-plane = 90°), whereas, the surface shear rate 
increases over the same domain. These two competing mechanisms 
appear to nullify each other around Q = 90 plane, thereby leading 
to a plateau In shear stress profile. The flattening of the shear 
stress profile becomes Increasingly evident as the value of n 
drops below 0.6. Furthermore, It is clearly seen In Figures 6.6a 
and 6.6b that for the range 0.6 < n < 1, the maximum shear stress 
occurs at & = 90^ plane, whereas, for highly shearthinning 
behaviour (n 2 * 0.4), the shear stress exhibits two maxima between 
70° < B < 110°, thereby giving rise to a cusp. This feature is 
clearly shown In Figure 6.6b on an expanded scale. 


Qualitatively 
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one can ascribe the formation of cusp to the aforementioned two 
competing mechanisms, i.e., the different rates of variation of 
viscosity and shear rate over the particle surface. A 

quantitative explanation for this phenomenon was sought by using 
the closed form stream f unct i on solution due to Kawase and 
Ulbrecht (1981a). However, a detailed analysis revealed only one 
extremum in shear stress at S = 90° plane It is not at all 

surprising as the analysis of Kawase and U lbr echt is applicable 
only for n close to 1 whereas, the cusp is observed for highly 
shear thinning behaviour, i.e., n % 0.4 No more explanation can 
be given for this behaviour at this stage. It is also of interest 
to note that as the vlaue of n drops further from unity, the 
maximum surface shear stress decreases. 

The surface shear stress distribution for a sparse assemblage 
(e * 0.9) Is shown in figure 6.7 for various values of n. The 
variation is seen to be qualitatively similar to that for « = 0 4 . 
The maximum shear stress decreases more or less uniformly with 
decreasing value of n, which can be explained by the gradual 
reduction In average shear rate on account of the large 
inter-parti cl e distance. 

Finally, the surface shear stress profile for an extremely 

dilute assemblage (in the limit of a single sphere) for « = 0.9999 

(R /R sg 20) is plotted in Figure 6.8. The maximum shear stress 
c p 

shows interesting trends with respect to the power law index (n). 

The value of t seesis to increase with power law index upto 
max 

about n & 0.8 Any further reduction in the value of n brings 
about a slow decrease in maximum shear stress initial ly , 
followed by a steep drop in r m&x * or ea n * Uni ike in 
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various 


th< ca.Be of dense assemblages , it appears that in the case of 

single sphere, the rate of increase of shear rate outweighs the 

corresponding decrease in viscosity, thereby resulting in net 

increase in the value of surface shear stress m the range (1 i n 

£08). As the value of n drops from 0 8 to 0.6, the two effects 

compensate each other, exhibiting near constancy for the surface 

shear stress. The only noticeable change seems to be a slight 

drop (less than 4%) in the value of t with a concomitant 

max 

increase in the value of shear stress near the axis This is in 
line with the discussion earlier made regarding flattening of 
shear stress profile with decrease in power law index. For lower 
value of n (< 0.6), the reduction in viscosity is much larger than 
the increase in shear rate. Therefore, r decreases rapidly 

with decrease in n . 

A further insight into the variation of surface shear stress 
under a range of conditions can be obtained by examining the 

variation of dimensionless surface viscosity and shear rate. 

Figures 6.9 to 6.11 show these results for a range of conditions. 
Evidently, as one would expect, the surface viscosity attains a 
maximum value on the axis (B = 0° and 180 ) due to the occurrence 
of lowest velocity gradients at these locations. Similarly, it 

assumes a minimum value at & - 90°. Recalling that in the case of 

Newtonian fluid, the dimensionless viscosity at the sphere surface 
is constant, it appears that in multi-part icle systems a major 
part of the particle surface is exposed to a constant value of 
viscosity (but much lower in magnitude than the Newtonian case). 
This prompts one to conclude that the variable viscosity effects 
are confined mainly to a small region in the vicinity of the axis 
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(& = 0° and 180° ) . 

For Newtonian and moderately shearthinning fluids, the 

surface shear rat e conforms to sinusoidal variation; departures 

are, however, observed for dense assemblages and highly 

shearthinning conditions (Figures 6.10 and 6 11). The surface 

shear rate also shows a cusp similar to that seen in shear stress 

profile. It appears that the interactions between the dependences 

of shear rate and local viscosity on n exert an appreciable 

influence m determining the value of t around & — 90° plane. 

max 

For instance, when n = 0.8 and « = 0.4, the drop m the value of 
viscosity from & = 0° to & = 90° is of the order of 60%, whereas, 
the corresponding change for « = 0.9999 is only 30%. This is 
precisely why the surface shear stress shows different types of 
dependence on n for dense and dilute assemblage. 

The trends of shear stress along with the surface shear rate 
and viscosity variations are helpful in shedding light upon the 
behaviour of friction drag coefficient (C^p) ior various voidages 
and power law indices. For instance, at low voidage the friction 
drag decreases monotonical ly with sharp drop near n = 1.0 and 

slower drop at lower values of n. At a higher voidage (« = 0.9), 

the friction drag still decreases monotically, but at a slower 
rate. However, while approaching the single sphere limit (« = 

0.9999), C n „ manifests a mild increase In the beginning followed 
by a sharp decrease as n is progressively reduced. These trends 
are in exact correspondence with the overall behaviour of shear 
stress profiles discussed earlier. 

The pressure drag coefficient is likewise determined by the 
pressure distribution on the surface of the particle. Typical 
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surface profiles are shown in Figures 6.12 to 6.14. Note that in 
these figures, the value of pressure at 6 = 180° has been used as 
the reference point. For dense assemblages (in the packed bed 
range), the maximum surface pressure (& = 0°), is seen to decrease 
rapidly with the decreasing value of n, akin to the surface shear 
stress. Interestingly, for any two values of n, the ratios of 
maximum surface shear stress and surface pressure are nearly the 
same. It comes as no surprise since pressure and viscous forces 
balance each other in the absence of inertia. For the creeping 
flow of Newtonian fluids, the surface pressure is known to vary as 
cos & ; however, a departure from this behaviour is observed for 
highly shearthinning conditions, showing a milder pressure 
gradient in angular direction than suggested by cosine dependence 
m the neighbourhood of B = 90° plane. As the concentration of 
assemblage is reduced (<e = 0.9), the variation of pressure with & 
for different values of n is more or less uniform (See Figure 
6.13). Here also, the behaviour is at variance from cosine 
dependence for low values of n. 

6.1.3 Flow Field in the Ceil 

Typical dimensionless pressure and velocity profiles within 
the cell are shown in Figures 6.15 to 6.20 for various values of n 
and «e An inspection of figure 6.15 suggests that in the packed 
bed region (e = 0.4), the pressure decreases primarily due to the 
variation in & direction for all values of n considered. In the 
radial direction, the pressure variation is significant only in 
the region adjacent to the axis , where the flow is predominantly 
radial At 6 = 90°, there is no radial variation of pressure at 
ny voidage or power law index, since the flow i© purely in & 
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Fig. 6.16 


Non-dimensional pressure distribution m the cell 
a t g = 0.9. Effect of shearthinning (n = 0.8 and 

0.4) . 



Non-d imensional radial velocity. 



Fig. 6.17 Variation of non-dimensional radial velocity 

(V /U ) in radial and angular directions at 
r * co 

e=0.4 and n = 


0 . 8 . 
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direction at the raid-plane The small radial variation of 
pressure at low voidage is a consequence of email available radial 
gap and this feature has similarities to the usual approximations 
invoked in boundary layer theory. It is seen that the rate of 
change of pressure in radial direction is larger near the particle 
surface than at free boundary of the cell for all angles. This is 
simply due to the fact that as one traverses from the particle 
surface to the cell boundary, there is gradual reduction in radial 
flow which eventually vanishes at the cell boundary and hence 


(tfp/^r) goes to zero at the cell surface. 

On the other hand, in case of dilute assemblages (« = 0.9), 
the variation of pressure in radial and angular directions is of 
comparable magnitude (See Figure «.1«) This is due to the fact 
that distance travelled by a fluid element in these two directions 
is of the same order. At. « - 0 . 9 also , the drop m pressure in 
radial direction is higher near the particle surface than the drop 
near the cell surface. For Newtonian and mildly shearthinning 


fluid behaviour, the pressure varies monotonlcal ly over the entire 


a. hAr . e seems to be an adverse pressure 
flow domain, whereas , there seem 

.. + SU rface for highly shearthinning fluids, 
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Fife 6.19 Non-dimensional angular velocity (V^/U^) variation 
in radial and angular directions at « = 0.4 and 


n = 0.8. 
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Fig. 6.20 Non-dinenaional angular velocity C W variation 
in radial and angular dlrectiona at « 0.9 and 
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symmetrical about & = 90 plane which is typical of creeping flow 
The elope of radial velocity curve is zero on the particle 
surface. Using cont inuity equation in conjunction with the 
boundary conditions (V r = cos 6 and V & = - sin 0), it can readily 
be shown that SV^/Sr is indeed zero at the particle surface. 
Consequently, the normal stress consists of only pressure and the 
friction drag arises only due to the shear stress. It is also of 
Interest to note that the radial velocity conforms to a cosine 
dependence on angle at any radial position within the cell. 
Qualitatively similar trends are observed for other values of 
power law index (1 £ n £ 0.2). 

Typical angular velocity profiles are displayed in Figures 
6.19 and 6.20. From the particle surface towards the cell 
boundary, V undergoes a change of sign from -ve to +ve somewhere 
close to the particle surface at all angular positions. Also, 
complete fore and aft symmetry is displayed. This feature would 
correspond to the streamlines originating from the particle 
surface, turning back somewhere near the particle and finally 
attaching themselves back to the particle in the rear. The radial 
location where V becomes zero would obviously coincide with the 

ff 

turning point of streamlines. At low values of voidage, the locus 
of such turning points of streamlines appears to be a spherical 
surface however , this is so only upto about n ~ 0.8 and with 

further drop in the value of n, it gets distorted into an 
ellipsoid. Indeed, some of these features can be clearly seen in 
streamline patterns depict ed in Figures 6.21 and 6.22. In 
general , the magnitude of at the cell surface is seen to be 
much higher than the part i c 1 e velocity except near the axis . Once 
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again using the overall continuity, it cau be demonstrated that V 0 
at the cell surface increases as ( 1 / 6 ) Further detailed 

examination of the results suggests that the velocity profiles are 
relatively insensitive to the power law index in dense 
assemblages. For dilute assemblages (« = 0.9), the angular 

velocity conforms to the expected sinusoidal variation only for 
Newtonian and mildly shearthinning fluids. For highly 

shearthinning conditions (n 0.4) and large values of voidage, 
the dividing surface takes an ellipsoidal shape. Another 

interesting feature exhibited by V. profiles is that at cell 
surface the velocity value remains more or less constant over a 
range of angular positions from 6 nc 60° to 120°. A reference was 
made in this regard while discussing the pressure profiles to the 
existence of adverse pressure gradient at the cell surface for the 
range of angles mentioned above. These pressure gradients seem to 
be partly responsible for V remaining constant around the mid 
plane. For all conditions of « and n, the radial gradient of 
angular velocity is the steepest at the particle surface and it 
decreases gradually as the cell boundary is approached. This is 
simply due to the condition of zero shear stress at the cell 
surface . 

©.1.4 Streamline! Patterns 

Owing to the difficulty of plotting, the streamline patterns 
are shown only for a dilute assemblage (« = 0.9) for two values of 
power law index in Figures 6 21 and 6.22. As noted earlier, the 
streamlines are made up of closed loops originating in the front 
and terminating in the rear of the particle. In addition, one can 
also discern closed streamlines forming a vortex adjacent to the 
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particle However, all streamlines, display complete fore and aft 
symmetry for both the values of power law indices, as expected in 
creeping flow regime Furthermore, the \ortex diminishes in size 
with increasing degree of shearthinning It is to be recognised 
that the vortex is sustained by the shear stress caused due to 
flow of fluid displaced by the moving particle. At lower values 
of n, the shear stress exerted by the displaced fluid is smaller 
due to the decrease in effective fluid viscosity. This, in turn, 
results in reduction of the strength and the size of vortex. It 
is also seen that as the value of n drops from 0.8 to 0.4, the 
centre of vortex shifts closer to the particle surface, thereby 
providing more area for fluid flow to occur. Concurrently, the 
tangential velocity decreases as has been noted in the angular 
velocity profiles (Figures 6.19 and 6.20). 


6.1.5 Comparison with Available Results 

It is of interest to compare the present predictions of the 
total drag coefficient with the other contemporary studies in this 
field. It is customary to express drag coefficients of creeping 
flow through multiparticle assemblages in the form of deviation 
from the Stokes drag. This is accomplished by introducing a drag 
correction factor (Y) defined as 


Y = 


C D Re 
2T” 


(6.5) 


The quantity Y is expected to be a function of n and « for 
multi -particle assemblages. The variation of Y with n and « is 
shown in Figure 6.23. 

Mohan and Raghuraman (1976a) obtained upper and lower bounds 
on the drag coefficient of assemblages for 1 £ n £ 0.5 and 
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Table 6 2 Comparison of the present numerical results with the 

available results on drag correction factor for 
creeping power law fluid flow. 
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o 99 > «= i 0.4 A comparison between their results and those of 
the present study, is shown in Table 6 2 An examination of this 
table shows that two values are in excellent agreement upto abhout 
n = 0.8 and as the value of n drops further, the two values differ 
increasingly. Kawaee and Ulbrecht (1981a) also presented a closed 
form expression for Y and their values are also included in Table 
6 2. Their analysis seems to be in reasonable agreement with the 
present study only upto n *£ 0 8 . It is believed that the results 
of the present study should be more accurate since they are 
derived from a full scale numerical solution of the equations of 
motion and continuity. On the other hand, the analyses of the 
previous authors were based on approximate solution techniques. 

Admittedly, numerous experimental studies concerning the flow 
of power law fluids through packed beds have been reported in 
literature but not all authors have provided the complete 
Information for their results to be recalculated in the form 
required here. Representative experimental results of some 
investigators were analysed and recast in a form which facilitates 
comparison with the present results as shown in Figure 6.24. All 
these results are pertaining to the packed beds with the value of 
voidage close to 0.4. The power law index of the fluid in these 
experimental results ranges from nearly 0.8 to 0.4. Each of these 
experimental results represents the average value of several runs. 
Keeping in view the idealizations involved In the present model 
such as the absence of wall effects, uniformity of particle 
distribution, etc and usual experimental errors In measurement of 
non-Newtonian flow through packed beds, the agreement between the 
experimental and present numerical results is seen to be 
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Fig . 6.24 Comparison of experimental results of other 

investigators with the present numerical results for 
flow through the packed beds. 
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satis factory 

Development of a general fluidization correlation and 
comparison with the available results on bed expansion 
characteristics of fluidized beds has been discussed m Sec. 6. 3. 

6.2 FLOW AT INTERMEDIATE REYNOLDS NUMBER 

The objective of the study reported in this section has been 
to identify the effects of inertia on transport phenomena 

associated with non-Newtonian fluid flow through multi-particle 
systems. This ef f ect i vely boils down to the analysis of drag 
behaviour and detailed flow field structure at several Reynolds 
numbers beyond creeping flow regime, with voidage and power law 
index as parameters. Although, numerical results have been 
obtained in this study covering a range of parameters (IS Re < 20, 
0 3 s: « < 0.9, 0.4 S n i 1.0), only drag results have been 
presented in all the cases investigated, for the sake of brevity. 
Furthermore, surface pressure and surface shear stress 
distributions have been analysed in detail for Re = 20 only. This 
value of Reynolds number was chosen based on the consideration 
that the inertial effects may begin to manifest in a prominent 
manner around Re = 20; also, for higher values of Reynolds number 
the free surface cell model is likely to become inapplicable as 
discussed in Chapter V. An insight into the nature of flow field 
prevailing around a particle in an assemblage is provided by the 
streaml ine patterns shown for two values of power law index (n— 0.8 
and 0.4). Finally, detailed comparisons with appropriate 

experimental results on flow past packed beds and single spheres 


are presented. 
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©•2.1 Surface Shear StrtcB and Pressure Profiles 

Figures 6.25 and 6.26 show the typical shear stress profiles 
on particle surface for two values of voldage (0.4 and 0.9) and 
for a range of values of power law index (1 * n £ 0.2) at Re = 20. 
Unlike in the case of creeping flow regime, the shear stress is 
seen to be asymmetric about the mid-plane = 90°) and attains a 
maximum value in the rear of the particle. It is reasonable to 
expect that due to Inertial effects in the intermediate Reynolds 
number regime, a mass of fluid will be moving along with the 
particle in the front region with more or less same velocity as 
the particle. This phenomenon is more clearly seen in the 
streamline plots presented in Figures 6.40 and 6.41. Such a 
feature is akin to the presence of stagnation region in the case 
of fluid flow over a fixed particle The stagnation region for a 
fixed particle occurs in the rear portion, causing maximum shear 
stress to shift to the front region. In the present case, 
however, the stagnant region occurs in the front of the particle 
and the maximum shear stress shifts towards the rear region This 
is because in Happel cell model, particles are assumed to be in 
motion while the fluid at large Is stationary. 

Uith regard to the influence of power law index(n), 
examination of Figures 6.25 and 6.26 reveals a systematic drop in 
the value of shear stress at all angular locations with decrease 
in the value of n. This is a direct consequence of the 
shearthinning behaviour. The variation in shear stress appears to 
have an approximate power f unct i on type dependence on n, showing 
sharp drop for mildly shearthinning fluids (n close to 1.0) 
compared to the strongly eh ear thinning fluids (n S 0.5). It is 
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also interesting to note that even at Re = H , the ah *ar stress 
for Newt or ia i fluid (n = 1.0) varies in a nearly einuaoidal manner 
in dense assemblages 0*04) This is similar to the variation 
observed in creeping flow regime and is believed to be a 
cons equence of the strong influence of viscous force in the 
confined geometry of the multi-pat tide system, even at a fairly 
high Reynolds number of 20. With a decrease in the value of n, 
the shear stress profile becomes flatter, indicating a trend 
similar to that of the creeping flow regime At lower values of n 
0.6), the shear stress profile is found to be flat even beyond 
the midplane (© = 90°) In the rear portion, it slowly rises to 
reach maximum value and then suddenly drops to zero near the rear 
axis. The flatness of the shear stress profile in the front 
region can be attributed to the mutually offseting trends of the 
increase in surface shear rate and the decrease in surface 
viscosity with respect to angular position. Furthermore, due to 


inertial 

effects of the 

fluid 

in 

front 

of the 

particle, the 

velocity 

gradients are 

smal 1 

in 

that 

region. 

As described 


earlier, this leads to an asymmetric shear profile with maximum 
shifted towards the rear of the particle. 

As the inter-particle separation Increases (by increasing the 
bed voidage), the magnitude of the shear stress drops. This is a 
direct consequence of the reduction in surface shear rate, due to 
lower velocity gradients in the gap between the particle and the 
cell surface. The flattening of the shear stress profile with 
reduction in n is also pronounced at « = 0 . 4 , for the same reason. 
Comparison of the shear stress profiles corresponding to both 
voidages (« = 0.4 and 0.9) for a fixed value of non-Newtonian 






% f ace pressure (p f\ 

Non-dimensional suriace f ^s' 

. „ ara i values of power law index 

distribution at several values 

e = 0.4 and Re = 20. 


Non-dimensional surface pre ssure , 
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fluid behaviour index (n = 0.4), ehova that at higher voidage (« = 
0.9) the flat region is confined to 0° < <? < 60 ° f whereas, at 
lower voidage (« = 0.4), it extends upto © = 100°. Another 
important difference between the results of two voidages is that 
in sparse assemblages even the Newtonian shear stress profile does 
not follow the sinusoidal variation It is also observed that the 
location of maximum shear stress with respect to the angular 
position remains almost invariant with the change in the value of 
power law index. The implication of this feature is that the flow 
field remains insensitive to the variation in local Reynolds 
numbers at higher voidage. 

The corresponding pressure profiles for the same values of 
the kinematic (Re % 20) and physical parameters (<e 0.4 and 0.9) 
are presented in Figures 6.27 and 6.28 over a range of power law 
indices (0.2 < n ^ 1.0). For densely packed systems (« = 0.4), 
the pressure profile for Newtonian flow approximates a cosine 
profile, eventhough Re = 20 is well beyond the creeping flow 
regime. This affirms the contention that viscous forces are 
important for densely packed assemblages even at Re = 20. 
Furthermore, It Is observed that pressure profile increasingly 
deviates from cosine dependence on angle as the degree of 
non-Newtonian behaviour becomes more and more pronounced. Indeed , 
for low values of n, the surface pressure is seen to attain a 
minimum value in the vicinity of & - 90 plane. These features 
are qualitatively similar to that of a potential flow pressure 
profile. Thus, the occurrence of minimum pressure around 6 * 90 , 
is a consequence of flow acceleration around the aid plane. 
Indeed by combining the fact that Is maximum at aid plane with 
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the Bernoulli’* equation, it can be shown that the pressure goes 
through a minimum value around the mid plane. Between & = 90° and 
& = 180 , the pressure increases due to flow decleration 

However, the pressure recovery is not complete due to viscous 
dissipation of the momentum Thus for any given value of n, the 
rear stagnation pressure is lower than the frontal stagnation 
value. It is also interesting to note that the change registered 
In the value of surface pressure with the change in power law 
index, decreases as n Is reduced. Effectively, this means that 
when n Is changed from 1 0 to 0 8, the decrease In pressure Is 
more compared to that which occurs for change In n from 0.6 to 
0 4. Qualitative explanation for such a behaviour can be taken 
from the dependence of viscosity on power lav index discussed 
earlier in context with the creeping flow situation (Figure 6.9). 
As viscosity becomes very small at low values of n, the scope for 
further reduction in It remains limited Therefore, only minor 
changes occur In the pressure profile, with reduction In n at 
lower values of n . In the case of a hypothetical fluid with n % 
0, the pressure profile will correspond to that of an Ideal fluid 
flow occurring within the cell and such pressure profile would 
exhibit fore and aft symmetry about the mid plane. Indeed the 
present results are seen to approach this limiting behaviour. 

At high voidage (€ * 0.9), the surface pressure profile 

monot oni cal ly decreases from the front stagnation value to the 
rear with few minor exceptions. The monotonic decrease Indicates 
that viscous effects are significant for all power law indices 
considered (0.2 < n < 0 8) even at Re = 20. Infact, the velocity 
within the cell gap varies as (1/d), where 6 is the radial gap 
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between the particle and th 4 i 

ana the cell surface (See Appendix- 1 1 ) . The 

inertial effects are . ,„2. 

are therefore ~ 0(1/6 ) while the viscous effects 

are 0(1/6 ) Thus, at high voidage and Re = 20, the inertial 

effects are not as important as in the case of low voidage. It is 

also worthwhile to note that the change In pressure with the power 

law index is much more gradual in this case than seen in Figure 

6.27. This is due to a reduction in the average shear rate 

arising from the larger radial gap within the cell. At low values 

of n (:£ 0.4), the pressure profile is seen to fluctuate around & - 

90 plane Though the precise reasons for this behaviour are not 

immediately obvious, this may be attributed to the onset of 

inertial domination in the flow for given Reynolds number and 

voidage . 


6s 2. 2 Drag 

Pressure, friction and total drag coefficients numerically 
computed over a range of voidage (0.3 < <= < 0 9999), power law 
index (0.2 < n £ 0.8) and Reynolds number (1 £ Re £ 20) are 
summarized in Table 6.3. Figures 6.29 to 6.37 show the variation 
of the pressure, friction and total drag coefficients with 
Reynolds number and voidage for three values of the power law 
index covering moderate to high shearthinning characteristics of 
the fluid. A preliminary examination of all these figures 
suggests that the results conform to Re = constant (which is a 
trend displayed by creeping flow solutions) almost upto Re ss 10-20 
for all values of n and moderate values of e. For the extremely 
dilute voidages such as e % 0.9999 which are closer to the single 
sphere limit, however, the drag curve veers away from C D Re = 
constant around Re ss 5-10 or so. The proportional ityy of the drag 
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Table 6 3- Summary of drag (pressure, friction and total) 
coefficients for power law fluid flow past an 
assemblage of particles at intermediate 
Reynolds numbers 
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coefficient to flow Reynolds number reinforces the idea that for 
multipart! c 1 e systems, the viscous force plays a significant role 
and some of the features of the creeping flow solution are 
exhibited even upto Reynolds number of 20 . Strictly speaking, the 
critical value of Re at which the behaviour ceases to conform to 
Cp Re = constant, appears to be dependent on the values of and 
n Lower the value of voldage, higher i' a the critical Reynolds 
number, whereas, lower the value of the power law index, lower is 
the critical Reynolds number This observation in respect to the 
effect of voldage Is in qualitative agreement with the findings of 
El-Kaissy and Homsy (1973). In the light of above discussions, it 
is safe to conclude that the viscous effects continue to dominate 
the flow almost upto Re ft 20 m multiparticle assemblages under 
commonly encountered shearthinning conditions and voidage. 

Further examination of these figures reveals that the total 


drag coefficient draws a larger contribution from the pressure 
drag coefficient in the case of dense assemblages, whereas, 
frictional drag contributes in a greater measure in dilute 
systems. This reversal appears to occur around «= ft 0.7 or so. It 
can be shown from an order of magnitude consideration that for a 
small radial gap (<5) between the particle and cell surfaces, 
average angular velocity = 0(1/6) 


inertial force 
shear force 


= 0(1 /6 j 

= 0(1 / 6 2n ) 


and, pressure force 0(1/6 ) 

=2 for inertia dominated flows (i . e. low n, and/or high 
fl IT* i JR ** 

_ v , ,_ + i f or viscous flow (high n and/or low Re). It then 

Re) and m - x 1 u 

low voidages, pressure drag contribution will be 


follows that for 
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larger than friction drag co 6 £fiH. ft+ » ^ , 

e et t icient and vie e- versa . 

In the limiting case of a „ 

ot a single sphere (« = 0.9999), the 

preeaure drag mo.fflcl.nt ahowa anomalou. dep6ndenc « on voida([c 

und.r highly ahaarthlnnlng conditiona (Figure 6.31, l.a. pm.a.uru 

drag increases with increase, , a m . , 

ae in voidage. This is believed to be a 

conaaquanma of the aignlflcant changea In the flow field atruoture 

aa alngle aphere limit la approaohed, for an extremely 

ahear thinning fluid. One can poaelbly explain thla by noting that 

the flow la uaually circulatory in a multiparticle aaaemblage 

(Flgurea 6.36 and 6.39) whereaa the etreamllnea are non- 

circulatory In the caae of a single aphere This change in flow 

structure may cause the increase in pressure drop with voidage for 

an extremely dilute assemblage. The numerical values of r , C nr 

and Cp indicate that a power function type dependence (as 

discussed in Sec. 6. 1.1) is applicable for n 2: 0.5 even at Re = 20. 

However, the drag coefficient (Cp) near the single sphere limit (« 

~ 0.9999) is a continuously decreasing curve. Uith reduction in 

n, this behaviour is different from that observed for creeping 

flow at the same voidage. This difference In drag behaviour 

between the creeping and intermediate Reynolds number flow 

regimes, is evidently contributed by Inertia. 


6. 2. 3 Flow Field 

Figures 6.38 and 6.39 depict the typical streamline patterns 
(Re = 20 and e = 0.9) for two values of power law flow index. The 
flow is seen to be asymmetric in the gap within the cell about the 
mid plane. This asymmetry is clearly Indicated by a forward shift 
in location of the vortex in the direction of particle motion. 
The streamlines are seen to originate from the front portion of 








216 


the particle, turn back at a certain distance, and reattach 
themselves in the rear of the particle it is interesting to note 
that after leaving the particle surface in the front portion, the 
streamlines remain parallel to the axis, upto a certain distance 
away from the particle. In this region the stream function at a 
point can be roughly calculated as y = — y 2 , where y is the 

distance measured from the axis This is a feature of uniform 
flow with unit velocity. Therefore, a portion of fluid in front 
of the particle seems to be moving with the same velocity as the 
particle due to Inertial effects. This uniform flow region 
adjacent to the particle extends almost upto midplane (0 s» 90°) 
In addition to the streamlines attached to the particle, one can 
also discern small vortex like streamlines adjacent to the 
particle. As power law index is reduced the circulatory vortex 
adjacent to the particle shrinks in size As the vortex is 

sustained by shear stress in the fluid, it is hardly surprising 
that at a lower power law index (hence lower viscosity), the size 
or the vortex is small. However, it is to be mentioned that 
compared to creeping flow situation, the reduction in size of 
vortex with decrease in power law index is less apparent at higher 
Reynolds number (Re =20) 

6. 2. A Comparison with Available Results 

In order to ascertain the reliability of the present 
theoretical predictions, it is desirable to compare them with 
available experimental results in literature on flow of power law 
fluid past multiparticle assemblages at intermediate Reynolds 

« It is to be admitted that such comparisons are 

number - However, 

difficult due to the paucity of results at intermediate Reynolds 
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In a comprehensive study, Kemblovski and Mertl (1974) 
proposd the following expressions for the flow of power law 
liquids through the packed beds («= ~ 0.39 - 0.4): 

f bk = C150/Re bk ) + 1 75 K i / (V 2 C (K 2 - l) 2 + K 2 ) 1/2 (6 6) 

where , 

K. = K Re’ 

1 c bk 

Se bk' C < 2 V C1--0 

1-n 

V = ( ( 9n + 3)/n n (150 S«)~^~/12 

S = d 2 « 3 /150 (l-«) 2 

The quantities /j and f were found to be functions of flow 

c c 

behaviour index n only Values of the parameters fj Q and % for a 

constant n, were determined employing the method of least squares. 

Using the fifth order polynomials given by authors to evaluate p c 

and £ , it is found that when 

c 

n = 0.8 

? = 0.1924 

N c 

p = 0.8273 

c 

and when n = 0 6 

£ =0 1579 

c 

u = 0.7128 

c 

Equation (6.6) has been claimed to be applicable for both 
shearthinning and dilatant materials. However, the experimental 
data used in obtaining the empirical constants are restricted to 1 
> n > o.6. Therefore, the comparison between the present results 
and those calculated from Eqn. (6.6) is restricted to 1 £ Re SJ 20 
and 0.6 < n < 1.0 only. A summary of the results is shown in 
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®xamination of this Table suggests that the present 
theor ical results are in reasonable agreement with the 
experimental results (assumed here to be represented by Eqn 6 6) 
only upto Re ~ 10 and n ~ 0 8 The degree of divergence increases 
rapidly with the increasing non-Newtonian behaviour as well as 
with increasing value of Reynolds numbers. 

Table 6.4 Comparison of Theoretical Predictions with Experimental 
Results on Packed Beds («= = 0.4). 



n = 0.8 


n = 

0.6 

Re 

f bk (present) 

f bk (expt > 

f bk ( ' present) 

f bk (6Xpt 5 

1 . 0 

40.48 

45.87 

16.04 

22.891 

5.0 

8.1 

10.25 

3.25 

6.056 

10 . 0 

4.06 

6.16 

1.63 

3.903 

20 . 0 

2 06 

3.95 

0.86 

2 826 


In a recent study, Srinivas and Chhabra (1990) demonstrated 
that the values of the pressure drop for the flow of power law 
liquids through packed beds predicted by various empirical 
expressions available in the literature differ from each other by 
as much as a factor of 2. Furthermore , Sheffield (1975) 
re-analysed most of the data on the non-Newtonian fluid flow 
through fixed beds available in the literature, and concluded that 
the errors of the order of 60-70% are not uncommon. Bearing in 
mind these factors, the correspondence between the two results in 


f bk (present) 


3 _ 3 

4 C d« 
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Table 6 2 La believed to be satisfactory upto about Re ~ 10 and 
n ~ 0 8 For moderately shearthinning fluids (n = 0 6) the 
agreement with experimental results may be taken as satisfactory 
only at Re = 1.0. 

In order to examine the performance of the free surface cell 
model at intermediate Reynolds numbers for extremely dilute 
assemblages (« = 0.9999), the present results are compared with 
the available experimental results for a single sphere falling 
under the influence of gravity. At the outset, it is important to 
recognise that even the condition of « = 0.9999 denotes the 
CVV rat *° ^out 21, and based on the experimental results 
obtained in cylindrical tubes, the wall correction to fall 
velocity for this situation is of the order of 85; , which in turn, 
would contribute to drag coefficient in varying amounts depending 
upon the value of n. Thus, such a comparison must be regarded as 
only qualitative Figures 6.40 and 6.41 show typical comparisons 
for n = 0.8 and 0.6 with the experimental results of Chhabra and 
Uhlherr (1980). The degree of correspondence is seen to be 
qualitatively satisfactory upto Re ~ 20 for n ~ 0.8 and upto Rc < 
2.0 for n = 0.6. In concluding this section, it is to be 
mentioned that flow past packed beds and single spheres represent 
the two extreme conditions of particle assemblages and the 
intermediate voidage range is the characteristic of fluidized 
beds. In the following section, the performance of free surface 
cell model in dealing with the fluidized beds has been discussed. 



6.40 Comparison between theory (c = 0.9999) and experiments 
for single spheres (n = 0.8) at intermediate Reynolds 

numbers . 
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6.3 FLUIDIZED BED DATA AND COMPARISON WITH EXPERIMENTAL RESULTS 

A general fluidization correlation which may be useful for 

design and analysis purpose has been prepared essentially by 

extending the approach described by Zenz and Othmer (1960) for 

Newtonian fluid, to a power law fluid. This consists replotting 

the variation of drag coefficient vs. Reynolds number data for a 

2-n 2 n 

fixed value of n, in terms of Re^against (Re 1/n /C D )^ 

with voidage as parameter Such plots have been prepared for two 
values of power law index n (0.8 and 0 . 6 ) , and are presented as 
Figures 6.42 and 6.43 A plot of this kind is useful as it 
results essentially in a plot of particle velocity against 
particle diameter for a given value of n, with physical properties 
as parameters For a fluidized bed 


C 


D 


and Re 


. gd 

_4 P 

3 U 2 

CD 


(Pp-P) 

p 


p u 


2-n 

oo 


K 



( 6 . 11 ) 


( 6 . 12 ) 


From these two equations, it can be shown that 


2-n 

2+n 


Re 


2 

2+n 


= d a- 
P 


(6.13) 


n 


and [ Re 1/n /cJ 2+n = ^ ft (6.14) 

where ot and ft are solely dependent upon physical properties of the 

particle and the fluid and is the superficial flow velocity 

through the bed. Defining U as the free settling velocity of 

, U . 

single spherical particle, the relative velocity ratio [ |p J 


any voidage («s) can 


be obtained simply by drawing a line parallel 





d-C. 



Fig. 


6.43 Fluidization curve for n = 


0 . 6 . 



n 

to the |Re /cj n axis corresponding to the particular particle 
diamet ei Thus, it Is possible to plot the relative velocity 
ratio over a range of voidages (0 5 to 0.9) for a g •. fluid 
behaviour index (n) and flow conditions (Re) Such a plot is 
shown in Figure 6.46 for Re = 5.75 and n = 0.6 the selection of 
the rather unusual value of Reynolds number (Re = 5.75) was made 
because experimental data collected from several runs were 
available (Srinivas 1990) for these conditions These 
experimental points are also shown in Figure 6 44. The agreement 
between the theoretical predictions and the experimental results 
are satisfactory only upto a limited range of voidage (« :£ 0.6). 
At higher voidage, the divergence between the two results 
Iner eases . 

Another correlation to predict the mass and momentum 
transport to non-Newtonian fluids flowing through the fluidized 
beds has been due to Kawase and Ulbrecht (1985) The fluidization 
velocity-voidage relationship has been expressed as 



12.5 


18 F 1 (n) 


(9 




n ) 


i fn 


(i-*) 

€ (l+2n) n 


(6.15) 


where 


F 1 (n) 



(6.16) 


The predictions of eqn. (6.15) for fluid behaviour index, n = 0.6 
have also been shown along with the present results at Re Q * 5.75 
in Figure 6.44. It is evident that present results are in good 
agreement with the predictions of equation (6.15) over a range of 
voidage varying from 0.4 to 0.65. However, the divergence between 
the two results begin to increase beyond - * 0.7. The exact 



Relative velocity 


T 




FREE SURFACE CELL MODEL 

— KAWASE AND ULBRECHT (1985) 
A SRINIVAS (1990) 


0.4 0.6 

Void fraction, 6 


44 Relative velocity ratio (U^/U©) variation with voldage 
for Re = 5.75 and n * 0.6. Comparison with 

theoretical predictions and experimental results on 
power law fluid flow through fluidised beds. 



226 , 


reason for this behaviour is not obvious at this stage. It is tc 

be borne in mind that while calculating the relative velocity 

ratio in the present study, the free settling velocity of a single 

sphere in infinite medium has been substituted by superficial 

velocity in an extremely dilute assemblage (« = 0 9999). 

Furthermore, it can ,be readily shown that equation (6.15) is not 

valid at higher voidage range, and, in the limiting case of « 4 1, 

it predicts infinite relative velocity ratio which is physically 

incorrect Even at « = 0 9 and n = 0 6, it predicts U /U = 2.22, 

oo o 

whereas this ratio should always be less than unity. It is also 
appropriate to mention here that equation (6 15) is strictly valid 
for low Reynolds number flow, though It has not been mentioned 
exclusively by the authors However, since creeping flow 

characteristics persist even upto Re ^ 10, the comparison of 

present results for Re Q = 5.75 with the predictions of above 

equation does not entail any significant error The predictions of 
relative velocity ratio over a range of voidages (0.3 sS « ^ 0.9) 

are shown in Table 6. 5 at three different values of Reynolds 
number (0.1 , 5.75, 20.0) for n = 0.6 Also included in the same 
Table are the predictions of Kawase and Ulbrecht (1985) and those 
of Brea et al . (1976) which is essentially valid for low Reynolds 
number flow. 

On the basis of these comparisons, the present results are 
found to be satisfactory over a range of voidage (0.4 S « £ 0.6). 

Further conclusions regarding the performance of the model can be 
made only through exhaustive comparisons with established 

experimental results. 

The Influence of fluid behaviour index (n) on relative 
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velocity ratio velocity voideee relationship v.. also studied. 

The results for two values of fluid behaviour index (n = 0 8 and 

0 6) are shown in Table 6 6 at Re = 5 75 as well as at Re = 20 0 

° o 

It is evident that the reduction in the value of n from 0.8 to 0 6 
does not make any significant impact on velocity-voidage data even 
upto Re o = 20. This observation is consistent with the findings of 
Tonini et al . (1981). 

Table 6.5- Relative velocity ratio (U /U ) for moderate 

CD O 

shearthinning (n = 0.6) fluid at different Reynolds 
numbers-Comparison with available results 


R e o = 0 1 Re Q - 5 75 Re^ = 20 Kawase & Ulbrecht Brea et al 


(1985) (1976) 


0.4 

0.014 

0.0166 

0 0286 

0 018 

0 0077 

0.5 

0.04 

0 051 

0 0652 

0.0514 

0 .62520 

0.7 

0.19 

0 2296 

0.2437 

0.2947 

0.1504 

0.9 

0.619 

0 5925 

0 7391 

- 

0 5715 

0.9999 

1.0 

1.0 

1.0 

- 

1.0 


Table 6.6: Effects 

velocity 

through 

of power law 
ratio (U^/U o ) 
fluidized beds 

Index on relative 

for non-Newtonian fluid flow 


Re = 

5.75 


Re 

= 20 


o 



o 


€ 







n=0 . 8 

■o 

o 

II 

C 


n=0 8 

n=0 6 

0 4 

0.0173 

0.0166 


0.029 

0.02869 

0.5 

0.0446 

0.051 


0.043 

0.0632 

0.7 

0.1986 

0.2296 


0.2521 

0.2437 

0.9 

0.545 

0.5925 


0.6173 

0.7391 

0.9999 

1.0 

1.0 


1.0 

1.0 


6. A CONCLUSIONS 

Creeping flow is characterised by the dynamic balance between 
pressure and viscous forces. Thus both pressure drag and friction 
drag indicate a similar type of dependence on the non-Newtonian 
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peculiar features ar e ^ 

not observed Further , in the limit of 

approach to a single . 

sphere, the shear stress increases initially 

and then drops rapidly with reduction in n . For no .._ Newtonlan 

systems, the pressure profile is also modified from the cosine 

dependence on angular positions which is displayed in the case of 

Newtonian fluid flow T n f a . t „ „ , , 

intact, as n is reduced the pressure 

profile also becomes flat, which is markedly exhibited for the low 

voidage situation For the single sphere limit (« = 0.9999), 

surface pressure increases with decrease in n The observed 

trends of the surface pressure and shear stress with respect to n 

for a single particle case explain the existence of maximum drag 

which has been observed by several researchers A plot of surface 

viscosity variation indicates that a high value of viscosity 

exists near the axis Indeed, the non-Newtonian viscosity is 

nearly a constant on the particle surface, except near the axis. 

At intermediate Reynolds number flow some of the creeping 

flow features persist even upto Re as 20 for dense assemblages 

For instance, each of the individual drag coefficients (pressure 

and friction drag) is inversely proportional to the Reynolds 

number at a fixed voidage and power law Index. Also pressure drag 

drops mere steeply with voidage than the viscous drag, as in the 

case of creeping flow. A roughly power function type of 

dependence on n is also exhibited. The Inertial effects on the 

drag coefficients are noticeable only for high voidages, high 

Reynolds number and low values of n The surface pressure and 


shear stress profiles are significantly altered by inertia in the 
intermediate Reynolds number regime. The shear stress is very 
much asymmetric about the mid plane. The pressure profile shows 
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decease In value in the front portion and pressure re-., very in 
the rear of the particle This feature ia predominantly aisplayed 
at low values of n indicating that pressure and inert al forces 
approximately balance each other in such a situation. The 
streamline patterns also display asymmetry in the flow direction 
about the mid-plane of the particle 

Comparison of present drag results for creeping flow with 
upper and lower bound solutions of Mohan and Raghuraman (1976a) 
and perturbation solution of Kaw&se and Ulbrecht (1981a) indicates 
excellent agreement upto n *b 0 8 ; for n < o . 8 the divergence 
between the present results and those of earlier investigators 
increases. It Is believed that present results should be more 
accurate as compared to the other approximate analyses. 
Comparison with available experimental results on creeping flow 
through packed beds also shows reasonable correspondence. 

There is a paucity of both theoret ical correlations and 
experimental data pertaining to the bed expansion characteristics 
of fluidized beds involving non-Newtonian fluid flow A general 
fluidization correlation which may be useful for design and 
analysis purposes has been prepared in the present study by 
extending the approach described byy Zenz and Othoer (1960) to the 
flow of power law liquid- The relative velocity ratio (U^/l^) 
evaluated from this correlation for a range of voldages (0.3 < « SJ 
0.9) have been compared with the theoretical predictions of Kavase 
and Ulbrecht (1985) and experimental data of Srinlvaa (1990). A 
reasonable agreement is observed between present results and 
theoretical as well as eperimental results of previous 
investigators over the range of voidage, 0.4 S « < 0.6. 
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The present result* „» 

on v ^-r iation of drag coefficient with 

power law in r the case of «s = D 9999 have been compared, 

with the theoretical and exper imental results available for single 

sphere falling under the inf ly enc e 0 £ gravity, at low as well as 

intermediate Reynolds number®. For creeping flow (Re * 0 001) the 

present drag predictions are found to be closest to the lower 

bound values of Cho and Hartnett (1983) The maximum drag is 

found to occur at n - 0.6 For intermediate Reynolds number, the 

present predictions have been compared with the experimental 

results due to Chhabra and Uhlherr (1980) for two values of fluid 

behaviour index, (n 0.8 and 0.6) The degree of correspondence 

is seen to be qualitatively satisfactory. 
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SCOPE FOR FUTURE WORK 

In this work, it has been demonstrated that the use of as 
highly idealized as the free surface cell model can yield useful 
information about the transfer processes occur mg m 
multi-particle systems. Further refinements can be made by 
incorporating particle size distribution as often encountered in 
real life applications. The gap between theory and experiments 
can also be narrowed down by assuming more realistic shapes of the 
surrounding cell, e.g., cylindrical or distorted, etc. 
Alternately, one can model the flow in multi-particle systems as 
that in conduits of regular but complex cross-section. 

In this study, only shearthinning effects have been 
considered whereas most real materials display viscoelastic 
behaviour. The importance of viscoelastic effects in this flow 
field can be ascertained by replacing the power law model by one 
of the viscoelastic constitutive equation, e.g. , Maxwell fluid 

No consideration has been given to wall effects in 
theoretical developments whereas experimental results always 
entail a degree, howsoever small, of wall effects. These manifest 
in two ways: firstly, the container wall does provide additional 
surface area which comes in contact with fluid; secondly, the bed 
voidage approaches unity at the wall The latter can be taken 
into account by using voidage as a function of spatial 


coordinates . 
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Available experimental results differ widely even under 
nominally identical conditions, A part of this deacrepancy must be 
attributed to the incomplete rheological characterisation of test 
fluids, which have been seldom checked for viscoelastic and/or 
time dependent effects. Thus, there is a need for further 

experimental work in this area using well-characterized systems 
e & purely viscous, purely elastic fluids, etc Furthermore, the 
importance of wall effects can be ascertained by performing 
experiments with a range of particle to tube diameter ratios. To 
date, no experimental measurements on flow patterns prevailing 
around a typical particle in an assemblage have been reported 
Such measurements undoubtedly will further our understanding of 
the physics of flow. Thus, considerable scope exists for 

experimentalists as well as theoreticians for exploring the rugged 
terrains of flow in multi-particle systems 
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APPENDIX I 


ANALYTICAL SOLUTIONS FOR CREEPING NEWTONIAN FLUID FLOW PAST AN 
ASSEMBLAGE OF PARTICLES - DRAG COEFFICIENTS AND SURFACE PRESSURE 

DISTRIBUTION 


The dimensionless form of the Navi er-Stokes equation for 
steady, viscous and axisymmetric flow of Newtonian fluids, m 
terms of stream function (y>) may be expressed in spherical polar 
coordinates (r , B , ^) as 


Re f &~y & 


v 2 
E y> 


2 1 Sr SB K 2 2 . 

v r sin & 


. _ d . E 2 y 

; 36 Sc K 2 2^ 

r sm B 


) J sm B = E 2 Ew 
(A-l) 


where , 


S“ J sin B S , 1 <9 

~ 2 2 SB K sm & SB 

Sr r 


) 


For the case of creeping motion of sphere or creeping flow past a 
sphere, Eqn. (A-l) reduces to 


e 2 e 2 v = 0 


( A- 2 ) 


The radial and tangential components of velocity vector are 
related to streeam function as 

1 ( A-3 ) 


~~i — : 7 
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r sin & 

A general solution of Eqn. 
Brenner, 1965). 


(A- 4 ) 


(A-2) is of the form (Happel and 




(A x r + A 2 r + A 3 r + 


) sin B 


(A-5) 
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Therefore , 


V r = ' 2 c °s 6 (A^ 2 + A + tl + ) 

r r 3 

V © * Sln ® (4A r 2 + 2A, + ^ ) 

i 2 r 3 J 

r 

The four constants (Aj . , X< ) appearing 

evaluated using the following four boundary 
cell model 


(A-6 ) 

(A-7) 

in these equations are 
conditions of Happel 


At r = 1 ; V f = cos 6 and = - sin & 

and 


(A-8) 
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r SO 
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(A-9) 

These boundary conditions in conjunction with equation (A-6) and 
(A-7) lead to following set of equations, 


, i* i i 4 V , i - 5 


1A 1 + 2A 2 * A 3 
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•• A 44 jj 
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c R J 
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ZA 1 " ZT~ ~ 71 + ~ 5 ~ 
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= 0 


(A-lOa) 

(A-lOb) 

(A-lOc) 

(A-lOd) 


In equations (A-lOc) and (A-lOd), R c refers to the non-dimensional 
cell radius. For a given voidage and hence known value of R c , the 
above set of equations can be solved simultaneously to obtain the 
constants Aj , . A^ . In present work, these were obtained using 
Matrix inversion subroutines. However, these can also be obtained 
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explicitly as a function of non-dimensional cell radius through 
the following expressions 
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-R ( 2R + 1) 
c c ' 
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V 3 
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1 5 R 2 -0.5 R f 0 
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(A-13) 
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= 

1.5 A x + 0.5 A 2 + 0 25 

(A-14 ) 

where , 




f l 
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(A-15 ) 

f 2 
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(2R* 4 4R 2 4 6R c 4 3) 

(A-16) 

f 3 

= 

(2 R c 4 1) (R c - l) 2 

(A-17) 


Knowing the value of A^ A 2> A 3 and A fl the radial velocity (V f ) 
and tangential velocity (V fi ) at any point within the flow field 

can be evaluated 


SURFACE PRESSURE DISTRIBUTION AND DRAG COEFFICIENTS 

The surface pressure may be found by Integrating the & 
component of the equation of motion from the front stagnation 
point (© = 0°) over the particle surface i.e., r - 1. For steady 
flow, the 6 components of equation of motion in non-dimensional 

form may be expressed as: 
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(A-18) 
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At particle surface r equals unity Therefore, 


3p _ 2 f ^ V 5 „ 2aV 8 . , 2 * V r 

** ^ l a . 2 Sr a * 2 *> 


Y e A * W e 

—27 + cot * *r 

sm & 


(A-19) 


Substituting for and from equations ( A-6 ) and (A-7 ) , it can 
be shown that 


Sp sin , « tl v 

35 = — rr (40A i * 4 V 


(A- 20) 


Integrating eqn. (A-20) with pressure at & = 0° as reference 

pressure, (p = 0 at & = 0), it is found that 


40A + 4A 

P = Re { 1 - C ° S ^ 


(A-21) 


C = J - 4P sm & cos & d& = (40A^ + 4A^) 


(A-22 ) 


Furthermore, It can be shown that 
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Sr r 


, 0V - 

" - *3^- I sin 
r && J 


& d& 


(A-23) 


Substituting V and V g 


V from equation (A-6) and (A-7), it is found 


that 


c df = 5ES C “ 5A 1 + A 3 ) 

Finally, tha total drag coefficient le given by 


(A-24) 


C D _ C DP + C DF 


( A-25 ) 
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APPENDIX II 

ORDER OF MAGNITUDE ANALYSIS FOR SURFACE SHEARR STRESS 


For the streamline pattern set up in the cell, the 
dimensionless volumetric flow rate q created by the motion of the 
particle Inside the cell is equal to 
2 

q = n r U = n (B~l) 

as r = 1 and U = 1 on the particle surface. 


Denoting the dimensionless radial gap between the particle and the 
cell boundary as 6, one can express the area for flow at any value 
of & as : 


= fi(l + 6 ) 2 - n = n (2 + <5 ) <5 (B-2) 

The average angular velocity (V^) within the gap is therefore 
given by 

q __ ft 1 f pi — N 

V & ~ ' ESTW) = *('2+? T 

For dense assemblages (e.g <5 = 0.4), etc ), <5 would be small and 
thus , 


V 6 ~ 0 ( 1 /*) 

This can be used to deduce the surface shear rate as, 



r *** 

. ^ 

, sv _ 

A 1 r] 

i 

e = 

av 

[ Sr 

r 

r S6 J 

Sr 

r=l 1 


Utilizing (B-4 ) , e &v ~ 0(l/<5 2 ) (B “ 6) 

For Generalised Newtonian fluids, the average dimensionless fluid 
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lacosity at the particle surface is given by 


n~l 


n = 

ence, 

Y) - 


[^) ■ fr] 


n- 1 
~T~ 


. 2 (n- 1 ) 


CB-7) 


(B-8) 


"The dimensionless average shear stress can now be estimated as 


2 n * V S 2 n 
T Se -K— n 


Re " Sr Re _2(n-l) x 2 

o o 


,n 


’Or , t - 


Re<5 


2 n 


CB-9) 


"•which allows C^p to be written as 


C DF ~ Re ( ,2 ] 


n 


CB-10) 


Therefore, the friction drag coefficient can be expressed in a 
£ orm 


'DF 


Ax 


n 


where x = — 5 - - 

& L 

It is interesting to note that the velocity gradient or shear 
rate for the Newtonian flow situation is 0(l/<5 2 ). Therefore, 

C ~ 1 (e ) n 
U DF ITi ^ av newt. 

this observation can form the basis for evaluating the friction 
drag coefficients for non-Newtonian fluid flow problem from the 

corresponding Newtonian problem at the same voidage. 



/ £12542 




